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We study the scale-dependence of halo bias in generic (non-local) primordial non-Gaussian (PNG) 
initial conditions of the type motivated by inflation, parametrized by an arbitrary quadratic kernel. 
We first show how to generate non-local PNG initial conditions with minimal overhead compared to 
local PNG models for a general class of primordial bispectra that can be written as linear combina- 
tions of separable templates. Wo run cosmological simulations for the local, and non-local equilateral 
and orthogonal models and present results on the scale-dependence of halo bias. We also derive a 
general formula for the Fourier-space bias using the peak-background split (PBS) in the context 
of the excursion set approach to halos and discuss the difference and similarities with the known 
corresponding result from local bias models. Our PBS bias formula generalizes previous results in 
the literature to include non-Markovian effects and non-universality of the mass function and are in 
better agreement with measurements in numerical simulations than previous results for a variety of 
halo masses, redshifts and halo definitions. We also derive for the first time quadratic bias results for 
arbitrary non-local PNG, and show that non-linear bias loops give small corrections at large-scales. 
The resulting well-behaved perturbation theory paves the way to constrain non-local PNG from 
measurements of the power spectrum and bispectrum in galaxy redshift surveys. 



I. INTRODUCTION 

By far the leading scenario for explaining the 
initial seeds for the formation of structure in the 
universe is inflation, but we still don't know what 
kind of field(s) are responsible for it, nor their La- 
grangian. For simplicity it is usually assumed that 
inflation is driven by a single scalar field in its vac- 
uum state, that slow rolls and has a canonical kinetic 
term. Under these circumstances the fluctuations 
are very close to Gaussian, i.e. their reduced bispec- 
trum is of the order of the tilt in the scalar spec- 
trum [1]. The tilt is known to be rather small [2]. 
In the language of Eq. (1) below, it results in a non- 
Gaussianity /nl ^ 0.05, which is unlikely to be ever 
detectable. 

However, as discussed below, a large class of theo- 
ries can produce a Universe much like ours, but with 
small deviations from Gaussianity that are poten- 
tially observable with future experiments. Produc- 
ing interesting constraints on the physics in the early 
Universe requires searches for such deviations. The 
situation is analogous to particle physics, where two- 



point functions describe freely propagating particles 
and higher-order correlations (whose measurement is 
the goal of particle colliders) describe the physical in- 
teractions. In cosmology, measuring the power spec- 
trum only gives information on free fields in the in- 
flationary background. To learn about interactions 
and thus distinguish between models, it is essential 
that we measure higher-order correlation fimctions. 
In short, searching for primordial non-Gaussianity 
(hereafter PNG) is akin to using accelerators to con- 
strain the interactions in the early universe. There- 
fore, detecting primordial non-Gaussianity is our 
best chance to learn about the physics of inflation. 
The prospects for the next decade are very promis- 
ing [3]. 

Deviations from the assumptions enumerated 
above lead to observable non-Gaussianities. The 
simplest type of non-Gaussianity arises by going be- 
yond single-field models. A second scalar field, usu- 
ally called the "curvaton" could be light during infla- 
tion; such a field would come to dominate the energy 
density of the universe after the end of inflation, and 
then produce effectively a second reheating [4]. In 
addition, fluctuations could be generated during the 



reheating period when the inflaton energy density is 
converted into standard model particles with a fluc- 
tuating decay rate [5]. In all these models, which 
are very closely related [6] , non-Gaussianities can be 
characterized by a primordial gravitational potential 
at sub-horizon scales that obeys [7] 

^x) = cf>{x) + hr.{[Hx)? -{<!>')) (1) 

where is a random Gaussian field, leading to a 
bispectrum of the so-called local form, 

B'i''^' = 2/nlPiP2 + eye. (2) 

where = P^{ki) is the power spectrum, 'eye' de- 
notes cyclic permutations of the fc^, and /nl could 
naturally be of order /nl « 5 — 30 [8]. Primor- 
dial non-Gaussianity at this level should be de- 
tectable through measurements of the bispectrum 
of the CMB and in galaxy redshift surveys. Current 
constraints on such local shape are — 10 < /nl'^' < 
74 (95% hmit) from the CMB bispectrum in the 
WMAP7 dataset [2], and -29 < < 70 (95% 

limit) from the scale-dependence of bias in the power 
spectrum of galaxies and quasars in the SDSS-II sur- 
vey [9]. This constraint relies on the fact that local 
PNG generates a 1/fc^ scale-dependent bias, as first 
derived in [10], which justifiably has generated a lot 
of interest given the ever increasing volume of up- 
coming galaxy surveys^. 

Another possible deviation from the simplest infla- 
tionary scenario, staying within single-fleld models, 
is to go beyond canonical kinetic terms; the primary 
examples in this case are DEI inflation [12] and k- 
inflation [13]. In these models the bispectrum shape 
takes predominantly a form that peaks at equilateral 
configurations, the so-called equilateral shape [14], 
and a small contribution of a different shape called 
the orthogonal shape [15]. For the convenience of nu- 
merical calculations, these shapes are approximated 
by templates which are a sum of factorizable contri- 
butions. The equilateral template bispectrum may 



^ Of historical note, a quarter century ago [11] first pointed 
out that PNG can lead to a scale-dependent bias growing 
at large scales, when they proposed a PNG model with a 
trispectrum that leads to a 1/fe bias, but their result was 
somehow forgotten. 



be written as, 

O/NL 

+(Pi^p|P3 + cyc.) (3) 

where it should be noted that there are a total of 3 
cyclic permutations for the first term, and 6 for the 
last term. Here Pj = P^{ki) = Ak~^ is the power 
spectrum of the curvature perturbation (Bardeen 
potential), proportional to the Newtonian potential 
on subhorizon scales. The orthogonal template bis- 
pectrum reads instead [15] 

1 ^ortho ^ _3(p^P2 + eye.) - 8(PiP2P3)i 
6/nl 

+3(Pi^p|P3+cyc.) (4) 

More precisely, this template is only a good ap- 
proximation to the orthogonal shape away from the 
squeezed limit, otherwise more complicated tem- 
plates should be used [15, 16]. This issue is relevant 
to the calculation of the large-scale bias, as the more 
accurate template does not lead to a scale-dependent 
bias at low-fc whereas the simpler template in Eq. (4) 
leads to a 1/fc bias. On the other hand, such a be- 
havior is interesting from a phenomenological point 
of view as it is in between scale independence and 
the 1/k^ of local PNG. Moreover, such a scale- 
dependent bias is realizable in some models of in- 
flation [17]. We will leave for future study more 
accurate templates and proceed with Eqs. (3-4) in 
this paper. Current 95% limits from the CMB bis- 
pectrum in WMAP7 are -214 < /^l"" < 266 and 
-410 < f^'l^° < 6, respectively [2]. Limits on these 
shapes from galaxy surveys are, —419 < /nl"'' < 625 
and -179 < f^'l'^° < 6, respectively [18]. 

While these are the most often discussed devi- 
ations from the simplest inflationary physics, two 
more of the assumptions enumerated above may be 
broken. Relaxing the assumption of initial vacuum 
states leads to the so-called folded shape [14, 19], 
which gives 

1 ^folded ^ (p^p^ ^ ) ^ 3(PiP2P3)i 

d/nl 

-(Pi^p|p3+cyc.) (5) 



2 



Note that this shape ean be eonstructed as a hn- 
ear combination of the equilateral and orthogonal 
shapes. Thus, from the point of view of constrain- 
ing it in data, it is enough to consider the previous 
two non-local shapes. Finally, relaxing slow-roll can 
lead to more complicated shapes, as in resonant non- 
Gaussianity [20], which cannot be easily written in 
terms of the other shapes, due to its non-factorizable 
form, see also 

For a discussion of most of these models from effec- 
tive field theory see [21]. Note that in addition, there 
could post-inflationary PNG contributions, such as 
preheating and phases transitions, which we ignore 
in this paper (see e.g. [22]). 

The current constraints on primordial non- 
Gaussianity from galaxy surveys come from the 
scale-dependence of the bias in the power spectrum 
of biased tracers such as galaxies and quasars, which 
are sensitive to the local form (Eq. 2) of the pri- 
mordial bispectrum [9, 10, 23-25]. Adding the 
galaxy bispectrum, not only enhances the constrain- 
ing power of surveys, but also helps constrain other 
shapes (e.g. equilateral and orthogonal) where no 
strong scale-dependence of bias is induced. This in 
fact makes galaxy surveys in principle as powerful if 
not more than the CMB in the search for primordial 
non-Gaussianity [26-30]. 

One of the main arguments for being cautiously 
optimistic about large-scale structure (LSS) con- 
straints is that while CMB gains more 2D modes 
by going to small scales, where other sources of non- 
Gaussianity become important, LSS gets more 3D 
m,odes by going to larger scales where the physics be- 
comes simpler. But even if these forecasts were not 
achieved in practice, it is important to obtain con- 
straints on PNG from large-scale structure (LSS) for 
several reasons: i) the kinds of surveys needed to ac- 
complish this are already planned, mainly to study 
cosmic acceleration, ii) if there is a scale-dependence 
in /nl it is essential to have a wide coverage of scales 
and CMB+LSS does this very well, iii) something as 
fundamental as claiming that the initial conditions 
are non-Gaussian will require confirmation by an in- 
dependent method. 

In this paper we want to implement N-body sim- 
ulations initial conditions for the different models 
discussed above and introduce an efficient algorithm 
for doing so. In addition, we provide more rigor- 
ous theoretical predictions for these models from the 



peak-background split and contrast these with mea- 
surements in numerical simulations, previous PBS 
predictions in the literature and those of local bias 
models. 

A guide to the main results may be useful here. 
For readers interested in the algorithm for generat- 
ing PNG initial conditions of non-local type (Sec- 
tion II), Eqs. (34-35) give our implentation for equi- 
lateral and orthogonal templates. Appendix A gen- 
eralizes these results beyond scale- invar iance and 
Appendix B gives a self-contained prescription for 
more general templates, including PNG. For 
those interested in the peak-background split main 
result on scale-dependent bias in /nl models, it is 
given by Eq. (76). Table I summarizes our nota- 
tion for bias parameters. Section III E discusses how 
this result differs from previous PBS calculations in 
the literature, and Section IV discusses comparison 
against numerical simulations (see Figs. 5-8). 

II. METHOD FOR GENERATING A 
DESIRED BISPECTRUM 

It is well known that the first term in Eqs. (3-4) 
is easily generated by writing the Bardeen poten- 
tial $ as a local nonlinear function of an auxiliary 
Gaussian potential </>, Eq. (1), and we would like to 
write similar transformations from ^ to $ for the full 
expressions in Eqs. (3-4). Since the bispeetra in all 
cases scales as total power squared, we are interested 
in a quadratic non-local kernel K such that, 

^ = cP + hLK[cj>,<P] (6) 
which in Fourier space can be written as 

$(fe)=0(fc) + J /nl [5u]K{ki,k2)cP{ki)^{k2) 
£kid^k2, (7) 

where the Fourier-space kernel K is dimensionless, 
and [5u] = fc(fe — ^12) (fci2 = fci -|- fe2)- For sim- 
plicity we assume that /nl is a constant indepen- 
dent of scale, although there are theoretical reasons 
to expect scale-dependence [31, 32]. It is however 
straightforward to include scale-dependent PNG, see 
Appendix B for the generalization. 

The bispectrum of such a model reads simply. 
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= 2hLKi2PiP2 + eye. (8) 

where K12 = K{ki,k2), which from Eq. (7) has 
the exchange symmetry K12 = K21, i.e. fei O ^2. 
The goal to generate initial conditions is to invert 
Eq. (8) to find the kernel from a given bispectrum. 
Such inversion is not uniqnc because the kernel de- 
scribes pairwise couplings between two modes, while 
the bispectrum is the result of the "average" (sum 
over cyclic permutations of two out of three modes 
in Eq. 8) over all such possible pairwise couplings. 
Therefore, there is no unique way of constructing 
the pairwise couplings themselves (K) using only the 
"average" information provided by the bispectrum. 

One solution to this dilemma is to declare that the 
pairwise couplings are equal to the average regard- 
less of the pair of momenta, that is, use the stan- 
dard reduced bispectrum, as the kernel, i.e. K12 — 
B / {P1P2 + P2P3 + P3P1) (see [33]). However, while 
this inversion provides a unique kernel for a given 
bispectrum, it introduces other issues. This pre- 
scription has the much more restrictive symmetry 
under cyclic permutations of ^1,^2,^3. This means 
that two nearly antiparallel high-frequency modes 
fci sa — A;2 couple with the same kernel than a low fre- 
quency kz = —ki — fe2 to a high-frequency mode (fci 
or fe2), which seems a rather strong assumption par- 
ticularly in the squeezed limit (where fci w ^2 ^ ^3). 
See Appendix E for more discussion on this for the 
kernels we obtain in this paper. 

Furthermore, using such "average" kernel will in 
general give an incorrect configuration dependence 
to the "snake" diagrams of the four-point function 
that depend on the square of the quadratic ker- 
nel (usually characterized by the tnl parameter, 
see [34]), since the averaging and squaring do not 
commute^. While one can of course fix this problem 



^ The "snake" trispcctrum amplitude is, from Eq. (7), 
Tsnake = (2/nl)2 [K(ki2,-ki)ft:(ki2,k3)PlPl2P3 + 

11 perm.], where ki2 = ki -|- k2 and Pij = P{kij). Us- 
ing the reduced bispectrum as a kernel gives a non-trivial 
constraint between this amplitude and the bispectrum and 
power spectrum. A well-known example of non-local (but 
non-primordial) NG is that generated by gravitational evo- 
lution from Gaussian initial conditions [35], in which this 
constraint is violated. It would be remarkable if inflation- 
ary PNG satisfied this constraint. 



by adding a suitable cubic term 0{f^-^L[(j),(p,(j)]) 
which through the "star" diagrams cancels the snake 
contribution (and adds the appropriate one), this is 
somewhat unnatural and would be costly numeri- 
cally (as fixing the snake topology through star di- 
agrams cannot be written in separable form). It 
seems more desirable to have freedom at the level 
of the quadratic kernel that might be used to repro- 
duce the correct bispectrum and the snake-topology 
contributions to the higher-point functions. Further- 
more, one might be able to do so while keeping the 
kernel a sum over separable contributions, and thus 
considerably speeding up the generation of initial 
conditions. Here we provide a first attempt (ignor- 
ing four- and higher-point information) along these 
lines. 

In this paper we start from the bispectrum tem- 
plates and work to invert Eq. (8). The main point of 
the bispectrum templates is to have expressions for 
the bispectrum which are factorizable functions of 
the ki, to speed up numerical calculations. It is sim- 
ple to construct an inversion procedure that keeps 
this basic property for the kernel, which leads to 
fast implementation of initial conditions generation 
for N-body simulations. Let's focus for definiteness 
on the terms that have fci -f-)- fc2 exchange symme- 
try. If one looks at the bispectrum contributions 
with this symmetry, there will be one solution for 
the kernel K12 that generates each such terms; each 
distinct solution corresponds to matching K12P1P2 
to the distinct bispectrum contributions of the form 
{P^Pi + P^Pf)P^~"~^ for distinct values of a, (3. 
This means that the number of solutions will equal 
the number of different exponents that the three 
power spectra are raised to, e.g. in the local model 
there are two exponents (one and zero) and thus 
there will be two separable solutions (see below). 

By performing linear combinations of the multiple 
solutions we construct a family of kernels that repro- 
duce to tree-level the desired bispectrum. We then 
impose regularity constraints from loop corrections 
to restrict the free parameters so that the tree-level 
results are not spoiled at large scales. While our ex- 
pression for the kernel in Fourier space is general, 
the expressions we present in the main text for the 
Bardeen potential $ in terms of the Gaussian field 
(j) are restricted to scale-invariant spectra. These ex- 
pressions are generalized beyond scale-invariance in 
Appendix A. 



4 



Since the templates we use (Eqs. 3-5) ean be writ- 
ten as linear combinations of local form and two ad- 
ditional terms, we will first discuss how to generate 
each of these terms separately. Depending on the 
type of PNG one is interested in, these generators 
can then be summed up according to the desired 
linear combination. 

Note that our methodology can be applied to any 
primordial bispectrum that can be written as sum of 
separable templates, i.e. gi{k\) 52(^2) 53(^3) with gi 
arbitrary functions, not ncc(^ssarily power laws, see 
Appendix B. For other approaches to generating 
non-Gaussian fields see [32, 36-40]. 



A. First Generator: Local non-Gaussianity 

Let's see how our method works in the simplest, 
well-known case of local NG. We require that 

2/NLifi2PiP2 + eye. = 2/NLP1P2 + eye. (9) 

Note that there are two ways in which the RHS is, 
say, symmetric under fci o k2 exchange, one is the 
term P1P2 (a = /? = 1), the other is the combination 
P3(Pi + P2) (with a = 1 and ^ = 0). Matching the 
contribution from K12 to the first leads to the trivial 
solution, K12 = 1, using the second one leads to the 
non-trivial one 



K 



12 



= ^(1- 
2 VPi 



+ 



(10) 



where _Pj = P12 = P{\ki + ^2]). This means that 
more generically one can generate the local bispec- 
trum using the kernel 



1 

P2 



(11) 



where u is a free parameter, to be constrained later 
from one-loop power spectrum considerations, as the 
resulting (tree-level) bispectrum is independent of u. 
For a scale-invariant spectrum Eq. (11) leads to a 
simple form for the Bardeen potential. 



$ = -Z- + /nl [(1 - «) 0' + w V-'a-i (0V'a<A)l , (12) 



where the operators introduced can be most easily 
understood in their Fourier representation, 

d(j>= y^^(f>{x) = j e-'^-^ k(t){k)d^k (13) 

and as usual, 

V-^A{x) = -J e-'*'-'^ (^) A{k) d^k (14) 

and the inverse operation 

d-^A = V-V-M = j e-'^-^{^A{k)d^k (15) 

Note that all these auxiliary fields can be obtained 
by efficient Fast Fourier transforms (FFTs). For all 
our generators, not just local, our algorithm only 
needs FFTs and basic arithmetic operations, avoid- 
ing cumbersome convolutions. 

The choice u = in Eq. (12) is often used to 
generate local PNG, and for a good reason, as we 
shall see soon. 

B. Second Generator 

For the second term in Eqs. (3-5), we have (choos- 
ing the arbitrary normalization appropriately) 

2/nl-F!:i2PiP2 + eye. = 6/nl(PiP2P3)^/^ (16) 

and in this case we are dealing with a single to- 
tally symmetric bispectrum (no choice regarding ex- 
change symmetry, i.e. we just have a = /3 = 2/3), 
thus there is a unique solution for the kernel. This 
can also be seen by rewriting Eq. (16) as 



Ki2 + eye. = 3, 



(17) 



where = KijiPiPj/P^.y/^ with P^,- = P(%), 

and thus since fcj are arbitrary, Kij = 1, which im- 
plies 



p2/3 
r^2nd _ ^3 
^12 



(PlP2)l/3 



(18) 
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It is then easy to see that the second term in 
Eqs. (3-4) is generated by the following nonlocal non- 
linear function for scale-invariant spectra, 



$ = </'-/NLV-2(a0)2, 

which leads to a bispectrum for $, 



(19) 



= hi^PiPl'^Pl'^ + eye. 



(27) 



independent of s and t. The results for all the kernels 
so far are general, while those for $(</>) assume scale 
invariance; the generalization of the latter is given 
in Appendix A. 



(20) 



D. Constraints from one-loop corrections 



C. Third Generator 

For the third term in Eqs. (3-5), we have 
2/NLi^i2-PiP2+cyc. = 2/NLPi/'^P2^^P3+cyc. (21) 

which contains three possible ways of arranging the 
fci k2 symmetry, namely a ~ 2/3, /3 = 1, a = 
1/3, /3 = 1, a = 1/3, /3 = 2/3, and thus there are 
three solutions, respectively 



K(2 



.1/3^^-1/3) 



Kii = PII\P-'/' 



+ P-/^) 



(22) 



(23) 



K{1^ = P3 (p-2/3p-l/3 ^ p-l/3p-2/3^ ^24) 

Therefore, there is a two-parameter family of ker- 
nels that generate the desired bispectrum, i.e. 



Kl^^ = {l-t-s)Ki^+t Kii + s Ki 



III 

2 



(25) 



In terms of the Bardeen potential, we have that the 
third term in Eqs. (3-4) is generated by the fol- 
lowing nonlocal nonlinear function (assuming scale- 
invariance) 



$ = + /nl (1 -t- 5)9-1(030) +f V" 



+sV-'a-i(V>50) 



which leads to a tree-level bispectrum for 



(26) 



While it is not surprising that a given bispectrum 
may be generated in different ways, leading to differ- 
ent higher-point functions at C(/nl)' O'^^ must re- 
quire that loop corrections do not spoil the tree-level 
results at large scales (e.g. the primordial power 
spectrum and the bispectrum) , and this may lead to 
constraints on the free parameters introduced above 
if the kernels are too singular in the infrared (IR). By 
one-loop correction, we mean specifically the correc- 
tion to the $ power spectrum due to the presence of 
PNG. And by tree-level, we mean the contribution to 
the $ power spectrum from the purely Gaussian part 
i.e. the cj) power spectrum. In order to preserve the 
large-scale behavior that one wants (typically some- 
thing close to scale-invariant), one should make sure 
the one-loop correction does not dominate over the 
tree-level one in the IR. 

It is easy to see that the requirement that the 
one-loop power spectrum does not alter the be- 
havior at large scales (i.e. we only accept correc- 
tions that diverge as or slower as fc — )■ 0) intro- 
duces non-trivial constraints. Consider for simplic- 
ity the local model. Using the kernel in Eq. (11), we 
obtain for the one-loop corrections to the Bardeen 
potential power spectrum, which at tree-level reads 



5P^ = 2{Ah^f X / (fq 



{l-uf u{l-u) 
q^\k — k^ 



\2k^J 



3\2 



k-q\^ + q^ + (JLY /^ (|fc-g| +g ) 



g3 |A;- 



(28) 



which gives a k ^ correction to the spectrum with 
a divergent amplitude, unless u = 0, that is, only 
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the "standard" local model kernel {K = 1) is al- 
lowed. In principle such singular behaviors can be 
cancelled by adding appropriate terms of 0(/nl'/'^) 
in the Bardeen potential, but we ignore this possi- 
bility here; as wc have no 0{f^-^) information, and 
unless required by symmetry this precise cancella- 
tion would introduce fine-tuning. 

For more complicated bispcctra with more than 
one type of term (e.g. as in equilateral and orthog- 
onal models) the constraints follow similarly (with 
longer expressions). The generic linear combina- 
tion of kernels that generate the desired bispec- 
tra given by Eqs. (3) and (4) to tree-level are (see 
Eqs. 11,18,25) 



For simplicity, instead of making this choice, we 
rather use u = which removes the operators pro- 
portional to V~^d~^ and simplifies our initial con- 
ditions generation. We will explore the impact of 
different choices of u and other possible templates 
on the four-point function elsewhere. Appendix E 
provides other motivations for the choice of u. 

E. Implementation 

Therefore the kernels we use are (assuming scale- 
invariance, their straightforward generalization be- 
yond this is given in Appendix A) 



EQ 



-3K 



loc 
12 



2K 



2nd 
12 



3K 



3rd 
12 ) 



KoKT = -9JriT - Si^iT + 9if- 



3rd 
12 • 



(29) 



(30) 



The one-loop power spectrum leads to the follow- 
ing results for these kernels: for the most dangerous 
term at low-fc we have an amplitude 



DC (w - 2s) (EQ and ORT) 



(31) 



and choosing s = u/2 automatically makes the 
amplitude vanish in both cases. For the fc~^ ampli- 
tude we have (after setting s = u/2), 



1 



oc 



(1 - 3tf (EQ) 
(4 - 9t)2 (ORT) 



and again these choices, t — 1/3 and t — 4/9, make 
the respective k^^ amplitudes vanish (another way 
to see this is to look at the squeezed limit of the ker- 
nels, see Appendix E). The only free parameter left 
after removing the k^^ and fc""* divergences is u. In 
principle this can be chosen to minimize (but not re- 
move) the amplitude of the one-loop correction, 
but there is no physical reason to do this as the tree- 
level result is more dominant in the infrared. The 
choice that would minimize such terms corresponds 
to 



1 



20 
21 



0.95 (EQ) 



u = 



|§ « 0.79 (ORT) 



$EQ = <^ + /nl [ - 302 ^ 4 a- 1 (0 50) 

+2V-2(0 VV) + 2V-'(a0)'l, (34) 



for the equilateral model and 



$ORT = + /nl [ - + 10 (0 90) 

-h8V-2(0 v20)-h8V-2(S0)2l, (35) 



for the orthogonal case. As in Eq. (1), one should 
subtract the expectation value of the quantity in 
square brackets to impose ($) = 0. As discussed in 
Appendix E this subtracts an infinite constant that 
is more IR divergent than in the local case (which is 
also infinite). In any case, the zero mode does not 

(32) enter into the physics that just depends on deriva- 
tives of However, requiring regularity for the ker- 
nel that enters into this expectation value may mo- 
tivate further constraints on the free parameters, see 
Appendix E for further discussion. 

Given these kernels wc; generate initial condi- 
tions for numerical simulations. In Appendix C we 
present tests of our algorithm that show that we re- 
cover the correct bispectrum for the Bardeen po- 
tential $ in the initial conditions for each of the 
operators appearing in Eqs. (34-35). Given the 
Bardeen potential, we generate initial conditions 
using second-order Lagrangian perturbation theory 
(2LPT), which provides a simple yet efficient way 
of minimizing transients from initial conditions [41] 

(33) that arise from using linearized Lagrangian per- 
turbation theory (Zel'dovich approximation). This 
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nal (/nl = —400) models run on a 2.4 h^^ Gpc box 
with 1280^ particles, using the Gadget2 code [43]. 
For each PNG model we thus have a total volume of 
166( Gpc)^, the largest to date, which will allow 
us to test theoretical predictions of large-scale bias 
to a greater accuracy than before. For such choices 
of /nl the skewness of the primordial density field 
is positive for the local and orthogonal case, while 
negative for the equilateral model. See [44] for a dis- 
cussion of higher-order moments in these simulations 
and mock galaxy catalogs built from them. 

In Figure 1 we show the difference in the matter 
density reduced bispectrum, 



{P1P2 + P2P3 + P3P1)' 



FIG. 1. The difference in dark matter reduced bis- 
pectrum Q from that in the Gaussian initial condi- 
tions case at redshift z = 0.97 for triangles with sides 
ki = 0.06 h Mpc~^ and k2 = 1.5 ki as a function of angle 
9 between fci and k2. The deviations seen in the N-body 
simulations agree with the expectations from linear per- 
turbation theory evolution of the primordial bispectrum 
(solid) for all three models considered: orthogonal (blue 
pentagons, /nl = —400), local (red triangles, /nl ~ 100) 
and equilateral (green squares, /nl = —400), from top 
to bottom at = 0, n. 



is important in particular when measuring higher- 
order statistics such as the bispectrum and the 
statistics of rare events, e.g. the mass function at 
high mass. For example, a 2LPT initial condition 
imposed at z = 49 (our choice for the simulations 
we present here) is equivalent when measuring the 
mass function at z = 1,0 to having run the same 
simulation with Zel'dovich initial conditions at an 
initial redshift of z = 2040, 4228 respectively [42]. 

Our simulations are part of the LasDamas (Large 
Suite of Dark Matter Simulations) collaboration [42] 
runs, extended to PNG models. The cosmological 
parameters are r2„j = 0.25, ilb — 0.04, 51a — 0.75, 
h = 0.7, Us = 1 and as = 0.8. For this paper we 
present results based on 12 realizations of local (with 
/nl = 100), equilateral (/nl = —400) and orthogo- 



in each of the three models from the Gaussian initial 
conditions case at z = 0.97 for triangles with sides 
fci = 0.06 ft. Mpc~^ and ^2 — 1.5 fci as a function 
of angle 9 between fci and fc2- The symbols (with 
error bars obtained from the scatter among 12 real- 
izations) denote the measurements, while the solid 
lines correspond to the predictions of linear pertur- 
bation theory evolution of the primordial bispectrum 
for each of the models (i.e. Eq. 36 with B and P 
given by their primordial values scaled by the lin- 
ear growth factor). We see a very good agreement, 
which is further evidence that the initial conditions 
in each case have been correctly generated. 

Our 2LPT-PNG initial conditions algorithm, be- 
ing a sum of separable terms for the kernel, is very 
efficient. For non-local models the initial condition 
generation takes only 35% longer than for local mod- 
els, which for iVpar = 1280'^ particles takes about 5 
minutes in 320 cpus. This is several orders of mag- 
nitude faster than summation over modes methods 
that use non-separable kernels recently proposed in 
the literature [38, 40], which scale as N^^^. (as op- 
posed to A^par In A^par comiug from the use of FFTs 
in our case). Therefore, we can generate non-local 
initial conditions with minimal overhead over local 
models for all Fourier modes in the simulation box, 
without being forced to restrict the number of modes 
with PNG to a low-fc subset as in [38, 40]. 
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III. THE PEAK-BACKGROUND SPLIT 



Excursion-Set Basics 



We now turn to a derivation of the expected clus- 
tering in generic PNG models. For this purpose, it 
is useful to briefly review the peak-background split 
(PBS) argument that allows us to calculate the bias 
of collapsed objects [45, 46]. We will comment on 
how our approach differs from other accounts in the 
literature below, and also contrast the PBS predic- 
tions with local bias models. One of the main goals 
of this paper is to compute the relationship between 
halo and matter overdensities, i.e. the bias factors 
(see Eq. 98 below for a precise definition). 

In the excursion-set formalism [47] , halo formation 
can be described as a random walk of the smoothed 
linear density field S as the smoothing radius goes 
from very large (infinitesimal variance cr^, and thus 
tiny S) to crossing the linear threshold for collapse 
6c at some finite smoothing radius (which defines 
the mass-scale of the object). The random walk 
follows from the fact that changing the smoothing 
radius or variance (which can be thought of as the 
time variable in the random walk) one allows into S 
smaller-scale Fourier modes which change 6 stochas- 
tically. The first (smallest cr^) crossing of the barrier 
described by 5c is equated with halo formation. 

Despite using a number of simplifying assump- 
tions (random walks about generic points rather 
than peaks in the initial density field, simplified dy- 
namics, etc; see e.g. [48] for a recent discussion on 
this) the excursion set framework is very useful as 
it provides a non-perturbative model (that includes 
e.g. exclusion effects at small scales [49]) and at 
large scales is known to reproduce, in the Gaussian 
initial conditions case, the local bias perturbative 
expansion with linear and nonlinear bias parame- 
ters [50, 51] which are in reasonably good agreement 
with numerical simulations [52-54]. 

Within the excursion-set formulation of collapsed 

objects (dark matter halos) of mass m, their mimber 
density per unit mass (obeying the standard normal- 
ization condition J {dn / dm) m dm = p) is given by 



( dn\ _ 
\dmJ 



= — dm Ilo{Ss,am,Sc)dds, 

J — oo 



m 



dat 



dm 



(37) 



where cr^ is the variance of the small-scale density 
field smoothed with some filter at spatial scale R 
(identified at first-crossing with the Lagrangian ra- 
dius of the halo) with m = Anpa^R^ /3. In Equa- 
tion (37) 

M^cc^m) = I ° no(5„(7^,<5e)d5« (38) 

is the probability of first crossing the linear-theory 
threshold for collapse 5c between "time" (t^„ and 
(7^ -|- da"^ and no(5s, cr^, (5c) is the probability dis- 
tribution for the small-scale density perturbation 
reaching 5g by "time" af-^^ given the initial condition 
that 5s = when cr^ = 0, and it satisfies the absorb- 
ing barrier boundary condition no((5c, cr^, (5c) = 0. 
For Gaussian initial conditions and a top-hat filter 
in A;-space no((5s, a^, 5c) = PG{.5s,a'^) - Pg{'25c - 
5s, (T^), where Pq denotes the Gaussian PDF. The 
subscript in IIq denotes that 11 is computed for the 
initial condition that 5s = when tjj^j = 0. 

Note that except to make contact with the lit- 
erature, we don't assume universality of the mass 
function, which requires 

no((5.,a^,^c) = i^( — ,— ), (39) 

in which case <t^Jx) is a function of = 5c /(Jm alone, 
in fact 



(40) 



where f{v) is the usual Gaussian factor in Press- 
Schechter theory for Gaussian initial conditions. 



B. Halo Bias and non-Markovian effects 

We are interested in calculating the clustering of 

objects (halos), which can be derived from their con- 
ditional mass function that describes how the abun- 
dance of halos responds to a large-scale perturbation 
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6£, i.e. conditional on the initial condition that at 
"time" (j| the density field was Sf . The Lagrangian 
number density perturbation of halos is given by 

^ ^ 9^ I-oo ^{Ss,(yL^c; Se, <jj)d5s ^^^^ 

where n((5s, cr^, Sc'.Si.af) is the conditional prob- 
ability for 6s with absorbing boundary conditions 
at 6c at time cr^ given that 6^ = 6^ at time a| 
when the smoothing scale was much larger (note 
that < cr"^). As we mentioned above, it is im- 
portant to stress that we use a zero subscript for 
unconditional quantities and we write 

Ilo{6s, 6c, al) = U{6s, 6c, a^; 0, 0). (42) 

Since we are interested in the large-scale limit of 
the bias, wc work in the regime where aj ^ o"^„ 
(scales much larger than the Lagrangian size of ha- 
los) and look for the relation between 6h and 6^. 
For reference, the Lagrangian size of halos of mass 
m = 10^^, 10^^ MQ/h is Rl = 7, 15 h-^ Mpc. In this 
case the calculation simplifies significantly because 
we can approximate, 

n(<5„ al, 6c; 6,, aj) « n(,5„ a^, <5,; 5t, 0). (43) 

Note that it is often used, or implicitly assumed, 
that one can write this as 

n(<5s, o-^, 6c; 6e, 0) « Uo{6s - 6e, a^, 6c - 6^), 

(44) 

where the last approximation only holds if non- 
Markovian effects can be neglected, i.e. the statis- 
tics of the random walk depend only on the total 
amount walked {6s — 6^) and to be walked {6c — 6s) 
at a given time so one can use a shifted barrier 
{6c 6c — 6e) from the unperturbed initial con- 
dition {6^ = 0). This is an important simplification, 
because it means the statistics of collapsed objects 
can be obtained for small large-scale 6^ by Taylor 
expanding the unconditional mass function (which is 
determined by Ho), and therefore the halo bias can 
be obtained from derivatives of the (unconditional) 
mass function. 



While Eq. (44) is routinely assumed and the bias 
computed from derivatives of the unconditional mass 
function, there may be corrections to the large-scale 
bias that will result from this. For peaks in a ran- 
dom Gaussian field, computing the bias from deriva- 
tives of the unconditional mass function does not 
hold in general [55], although it holds asymptoti- 
cally at large-scales. In [56] it is found that viola- 
tions of Eq. (44) for halos can alter the large-scale 
bias, but [57] find that the shift 6c - 6i in Eq. (44) 
should be replaced in the non-Markovian case by 
6c — ol6i where a depends on the filter and the power 
spectrum (a = 1 for a sharp- fc filter). This change, 
however, does not affect the Fourier space large- 
scale bias (as opposed to the counts-in-cell bias), 
which is predicted rather accurately by the unmod- 
ified Eq. (44) when compared to numerical random 
walks. In practice, though, there should be non- 
Markovian corrections other than filter effects (e.g. 
PNG itself, or physics that is not included). 

In simulations with Gaussian initial condi- 
tions, [53] tested how well the bias computed from 
the unconditional mass function agrees with the 
large-scale bias measured from simulations in real 
and Fourier space. In this case, however, the as- 
sumption of Markovianity (which converts the 6i to 
a 6c derivative) is not enough to do the numerical 
evaluation of the derivative and one needs to further 
assume universality of the mass function (to convert 
the 6c to a v = 6c/ (Jm derivative). Proceeding in 
this way [53] foimd deviations at the 5 — 10%, par- 
ticularly at high-mass. 

As we shall see, fortunately, we will be able to 
proceed quite far for /nl PNG without the need of 
a full description of non-Markovian effects, although 
the same is not true for jinl PNG. For further dis- 
cussion of non-Markovian effects see [47, 48, 56-60]. 



C. Modulation Bias from PNG in the PBS 

In the presence of primordial non-Gaussianity 
(PNG), different Fourier modes are no longer inde- 
pendent, and the probabilities 11 for the small scale 
density field fluctuations get modified by the pres- 
ence of a long- wavelength mode of the Gaussian field 
(j) as follows 



10 



U{5,, al, 6,; 5e, 0) ^ U[S,, a^cj)) , Cp{^) , S,; 0], 

(45) 

where Cp = {S'p)c are the cumulants of the smah- 
scale modes'^. That is, unhke the case of Gaussian 
initial conditions where only the mean becomes in- 
dependent, all higher-order cumulants now depend 
on the long-wavelength perturbation because small- 
scale modes are coupled to them through PNG. 

It is precisely this (/(-dependence that we derive in 
this subsection. For readers interested in skipping 
the details and going directly to section HID for 
the results on scale-dependent linear bias to leading 
order in /nl, the main results needed are Eq. (64), 
and Eqs. (69-70). 

To proceed, although we are primarily interested 
in /nl PNG, we will generalize the model to account 
for a cubic ^nl kernel, 

$ = <n + /NLif[</',<n]+5NLi[<n, (46) 

To avoid very long expressions, we will suppress the 
contributions of the (/nl kernel in most equations, 
except when the extension is a non-trivial general- 
ization of the /nl kernel (see Eqs. 56 and further). 
We will also take into account all orders in PNG 
parameters (/nlAnl) to arrive at very general ex- 
pressions (see Eq. 66 and below). 

If one is interested in first-order in /nl calcula- 
tions, computing the modulation of the small-scale 
variance a{(l)) by (p due to non-Gaussianity (linear 
in /nl) is all that is needed, while we can evaluate 
the third-order cumulant C3 (related to the skewness 
S3 by C3 — S3 CT'^) ignoring the modulation by the 
large-scale mode <j), since modulation of C3 by the 
large-scale ^ starts at the four-point function level, 
and thus is at least C(/nl) o^c C(5nl) for PNG mod- 
els with cubic kernels. Similarly, we can put C4 and 
higher-order cumulants to zero. However, for now 
we will keep effects up to 0{f^-^) or ©(gNL) to see 
how our formulas apply to more general combination 
of quadratic and cubic PNG. 



Here we assume that 11 depends on the cumulants of 5s, 
but [61] find sub-leading corrections that depend on deriva- 
tives of correlators, which we ignore for simplicity. 



To calculate the modulation effects, we perform 
the peak-background split (PBS) between small- 
scale ("peak", denoted by s subscripts) and long- 
wavelength ( "background" , denoted by i subscripts) 
perturbations. This split takes into account that 
the long-wavelength modes have been conditioned 
on to take values <?i^(fc) (that is, we only look at ran- 
dom walks that satisfy this large-scale constraint). 
The split is often done in real space, by choosing a 
top-hat filter at some scale R (larger than the La- 
grangian radius of the halo Rl)- This is the appro- 
priate constraint when dealing with co\ints-in-cells 
statistics (at radius R), but here we are interested 
in the power spectrum where all modes are weighted 
equally. Therefore, we need to impose a more nat- 
ural constraint for Fourier-space statistics, i.e. a 
sharp-fc filter. 

Another choice we must make is on what variable 
to do the split. In the excursion set framework, the 
random walk is performed by the full, non-Gaussian, 
density field, which is equivalent to constraining the 
non-Gaussian rather than (f>i{k). However, 

technically it is easier to do the split in (p since for 
a Gaussian field long and short modes are indepen- 
dent, and this is presumably the reason why all past 
work on the subject has done so. However, arguably 
the more physical variable is $. after all this is what 
is supplied to the equations of motion of gravita- 
tional evolution that describe the physics of halo 
clustering. We shall see that there arc two instances 
where constraining 0^ differs from but for sim- 
plicity we will do the split on (p where general results 
are easier to derive to all orders in PNG and discuss 
where appropriate how the results change when con- 
straints on $ are used instead. 

We thus proceed with a split on (p in Fourier space 

clij SC3;l6 /Cgplit ? 

,n(«)(q) = 0(q)e(q-fc,piit), 

<P^'HQ) = m^ikspnt-q), (47) 

where 6 is a step function, i.e. a sharp-fc filter. The 
splitting scale obviously must satisfy /Cgput^L ^ 1, 
where i?^ is the Lagrangian radius of the halo. When 
computing the power spectrum at scale k, we will 
choose fcgpiit = k, analogous to choosing the split 
scale to be R when computing counts-in-cells statis- 
tics at scale R. Note that the precise choice of the 
splitting scale does not appear in the calculation 
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until one computes one-loop corrections, e.g. one 
could have chosen e.g. fegput = '^k and get the same 
trcc-lcvcl results. Our choice means that when com- 
puting the power spectrum at scale k we have inte- 
grated out all smaller-scale fluctuations, and loops 
contribute a small renormalization of the tree-level 
PBS bias parameters that results in their running up 
to scale k (see Section V). The choice of a sharp- fc fil- 
ter rather than a top-hat split has the added benefit 
of avoiding WTuikR) factors that would otherwise 
appear in the power spectrum (with R poorly de- 
termined) when computing loops. Therefore, in our 
treatment the variance of the large-scale modes is 



am 



(48) 



The small-scale field will be smoothed with a top- 
hat filter at the Lagrangian scale of the halo, as this 
is the usual way of identifying collapsed objects, and 
will be introduced shortly. We are thus imagining 
random walks where the filter changes from sharp- 
k to top-hat as the Fourier modes added go from 
large to small-scale (the transition from one to the 
other filter should be done smoothly at ~ fc, but 
the details won't matter as long as kRi, <^ 1, which 
as we shall discuss in Section V is in any case the 
regime of validity of the perturbation theory for bias. 
We write for the small-scale density perturbations in 
Fourier space, 

6,{k) = M{k)^,{k) = M{k) [Uk) + fNL [</','/>]" 

(49) 

where the function M{k) relates the density fluctu- 
ations to the Bardeen potential $ 



M{k) 



2D{z)T{k) 



(50) 



where Jl^ and Hq denote the z = values of mat- 
ter density in units of critical and the Hubble con- 
stant, D(z) is the growth factor in the matter era, 
with normalization £) = a at early times, where a is 
the scale factor. T(fc) denotes the transfer function, 
and the gravitational potential obeying the Poisson 
equation at sub-horizon scales at redshift z reads 

$grav = -(^/a)T(fc)$. 

In order to properly perform the PBS at quadratic 
(and higher-order) we distinguish between the long 



and short wavelength limits of the kernel (which in 
general can be vciry different, see discussion in sec- 
tion II above, and Appendix E for examples). That 
is, in Eq. (49) the upper index in K denotes that 
the kernel is non-zero only when evaluated for large 
external momentum fc, i.e. specifically 

k'^\'I>A] = J SD{k-k,-k2) K^'Hk,,k2) 

X [(j)s{ki)(j)s{k2) + '2(f>s{ki)(l)e{k2)], 
d^ki d^k2 (51) 

where we have used symmetry in ki and k2 to write 
down two equivalent contributions from terms linear 
in (be, and K<^'^ is a high-pass filtered version of K to 
ensure that the only small-scale 0s(fei) and (f)s{k2) 
contributions included in the integral couple to a 
higli-fc mode (fc > fcgpiit), the missing contributions 
below the splitting scale are included in the large- 
scale density perturbation for which we have, 



Si{k) = M{k) Mk) + /nl kI'\<P, </>] 



(52) 



where the upper index in K denotes that the kernel 
is non-zero only when evaluated for small external 
momentum fc, 

Ki'\(p,<i^]= J 6D{k-k^-k2) K^'Hki,k2) 

X [(t)s{ki)(t)s{k2) + (t>i{ki)(t)e{k2)], 
d^ki d^k2 (53) 

where K^^^^ is a low-pass filtered version of K to 
ensure that the only small-scale modes (psiki) and 
(j>s{k2) that contribute in Eq. (53) are nearly op- 
posite (and thus couple to a large scale mode fc 
below the splitting scale). By definition we have 
K = K^^^ + K^^^ . We can make all this explicit by 
writing (see Eq. 47), 

i^W(fci, fc2) = K^kiM) e(fci2 - fc.piit), 

jrW(fei, fca) = K{kiM) e(fcspiit - fci2). (54) 

Note our split at quadratic order automatically in- 
corporates all contributions to large and small scale 
perturbations, ie. we do not assume e.g. for local 
PNG that 

V^$« « VV« + 2/nl(<A^VV^ + • V<t>s) 

« V2</., + 2/nl<^^VV« (55) 
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the latter often "justified" by saying that V(/)s ~ 
at peaks. However, since we are deahng with peaks 
in the density then Wcps is proportional to the ve- 
locity field, which does not vanish. Eq. (55) is an 
approximation at the level of M (which involves the 
Laplacian), and there are similar approximations of- 
ten made about the filter function that defines the 
mass scale of objects. Our expressions automatically 
keep all such terms, we shall see that they contribute 
to scale-independent bias for local PNG (see also Ap- 
pendix D for more details). 

Let us now comment on using (f> rather than $ (or 
S) as the split variable. One way in which differences 
can arise is when the kernel K^^^ is singular enough 
in the squeezed limit that there are non-negligible 
contributions to the large-scale modes (or 6() in 



Eq. (52) coming from (pg through the first term in 
Eq. (53). When this happens, constraining on ((>( and 
can differ, as cannot be well approximated as 
a function of just (f>i, that is, fixing (j>i is not fixing 
the large-scale density field because the small-scale 
modes 0s make a non- negligible contribution to it. 
We will find an example for the equilateral template 
where this can lead to different conclusions about 
the scale-dependent bias. 

Prom Eq. (49) we can now calculate all the needed 
modulations of the cumulants of the small-scale fluc- 
tuations, which will now be space-dependent due to 
the background of large-scale modes that breaks sta- 
tistical translation invariance. For the variance after 
introducing smoothing on mass-scale m, keeping up 
to 0(/2l) and 0(5nl), 



a^{(l>e)=al + j P^{p)d^p [4/nlM„(p) ^j,(a;) + 4/^^ \^p{x)\^ + 6 g^^ M^ip) ^^j^\x) 



(56) 



where we neglected a —aj contribution since we assume af-^ ^ (t|, but this can be taken into account if the 
observation scale k is not too far from the Lagrangian radius of the halo Rl- Although it is often argued 
that in doing the PBS the split scale must be much larger than Rl, this is not the case. If URl is not much 
less than unity what happens is that the modes that have been conditioned on contribute to the halo as 
well, but that's not a problem as long as these conditioned "long" modes are later averaged properly when 
computing halo clustering at k. What is difficult when kRL ^ 1 is to do this averaging perturbatively, see 
Section V. 

In Eq. (56) the first and second-order ipp fields are given by 

Vp{x) = J Mm{p - q)K^'\-p, q)M<l>-"'-''d^q, (57) 

ip^^^x) = J Mm{p-q^2)L^'H-P,Qi,Q2) x Mcii)MQ2)e-"'''-^d^qi(i'q2, (58) 

with L being the cubic kernel describing (/nl non-Gaussianity (see Eq. 46). In Eq. (56), 

al = J Mm{kfP^{k)d^k + 2/2^ j Mm{kf[K^'\k - q,q)]^P4,{\k - q\)P^{q)d^kd\ (59) 

is the variance of the small-scale density fluctuations including non-Gaussian corrections and we have included 
filtering into a redefinition of M, 

M„(fc) = M{k) WT^ikRL), (60) 

where Rl is the Lagrangian radius of halos of mass m and Wth is the Fourier transform of a top-hat window. 
Similarly, we can write for the third cumulant. 
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c(3)(<^^) = c(^) + J P^{ki)d^kiP4>{k2)d^k2 \lAfl,^M^{kr)Mm{\ki+k2\)K'^'\kiM)^kAx) 

+ 185NL M„(fci) M^{k2) Vk.kA^)] = + ^c^^\cl>i) (61) 

where is the third-order cumulant calculated to the desired order in PNG (as in Eq. 59 for cr^), and 

iPk^kA^) = j L^'\qMM)Mm{ki2 + q)cl)e{q)e-'^-^d''q, (62) 

is analogous to ^pk but for the ^nl kernel. All these </7-fields can be written in terms of functional derivatives 
of small-scale density modes with respect to (fig and admit a simple diagrammatic representations which help 
derive these expressions. 

What we need to compute the bias factors are the functional derivatives of these modulated cumulants 
with respect to large-scale fields in Fourier space, where we want the bias (as opposed to real space). For 
the variance we have, 

(^^)^=e-'^-=^4/NL j Mm{p)Mrr,{\p-k\)K^^\-p,k)P^{p)dh= j B^^{q,k-q-k)d\ 

(63) 

where (. . .)o means evaluating at 0^ = and wc introduced the cross-bispcctrum between small-scale 
smoothed density perturbations 6{k) = Mm,{k)^{k) and (p. This result should not be surprising, as it is 
expected from linear response, i.e. 

\VMk))o - P^k) ' y^^>^n-n - - j ^ss4> 

Note that because (/i(fe) is Gaussian, there is no contribution to the cross bispectrum that depends on K^^\ 
only K^**) appears, and the two contributions can be written identically after a change of variables, as in 
Eq. (63). This can be traced back to the fact that we conditioned on the (j)^ modes rather than the large-scale 
non-Gaussian (or 5t). This in contrast with local bias models, which predicts primordial non-Gaussianity 
enters only through the small-scale density bispectrum. We'll come back to this later. 

To make our expressions more compact in the generic PNG case (i.e. beyond quadratic /nl kernels) we 
will also use the second cumulant notation for the variance when convenient and slightly abuse notation for 
p = 1 to denote 5(,, i.e. 

c« = <5,, c(2) ^ a2(</.,), c^P)^{[5s{<t>e)Y')c, c^^^ = c(^^^(^?>c, (65) 

thus one should keep in mind that p = 1 is special (in the sense that derivatives with respect to it will be 
related to the conditional mass function in the expressions below, as opposed to unconditional). 

To compute linear bias to leading order in /nl, Eq. (63) is all that is needed. To go beyond this, we must 
compute higher-order derivatives of the variance or higher-order cumulants. Rather than do case by case, 
we quote the most general result {p,q > 1), 

( -D'c^"^ \ ^ {c^^'Kr) reiki) ■.■re{k,))c 

\V(t)e{ki) . . . V(Pe{kg)Jo P^{ki) . . . P4,{kg) 

^ e-ife,,,.x J T!f +Jl Jq^, . . . , fei;, - -fei, . . . , -fe,) d\i . . . d^g^.i 

~ P^(fcl)...P^(fcg) ' 

^e-^^i^'-'^ 7j,g(fei,...,fe„m) (66) 



{a\ct>^)re{k)) = e-^'^-^ / Bjj, (q, k-q, -k) d'q (64) 
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where ki-q = "Yjl^i^i^ ''■^d T*^") is the primordial n*^ polyspectra, i.e. T'^'^^ is the bispectrum, T'^^^ is the 
trispectrum. Note the subindices of T'^p+i) denote that this expression involves the polyspectrum between 
p 5- fields (at large momenta) and q (/(^-fields (at low-fc). For p = \ there is no smoothing, thus 5 is replaced 
by (5, according to Eq. (65). Note that in deriving this result, it is crucial that 0^ is a Gaussian field, we will 
come back to this in Section IIIF. 

For example, to calculate linear bias at 0(.9nl) one needs the modulation of the third-cumulant by one 
large-scale field, i.e. p = 3 and q = 1 in the above formula, which involves the primordial trispectrum. One 
can compute directly the derivative from Eq. (61) and check that the appropriate (with all combinatoric 
factors) trispectrum given by Eq. (66) is obtained. Similarly, to compute quadratic bias at ^^(/nl) C(,9nl)j 
one needs the modulation of the variance by two large-scale fields {p = q ~ 2 m. the above formula) , and this 
involves a different trispectrum, which can be checked directly by taking derivatives from Eq. (56). Lastly, 
we note that for p = 1, Eq. (66) is proportional to the PNG kernels at order q, 

In{k,m) = M{k), = 2M{kr2) K^^KkiM), hg = q\ M{k,,g) K^-^^ku ■ ■ ■ ,kg) (67) 

We are now ready to calculate the bias for arbitrary PNG. We introduce shorthand notation to distinguish 
between real-space (where the excursion-set expressions were written) and Fourier space (where we want to 
calculate the bias), otherwise confusion can arise. Let 

A = A{x), A = A{k), (68) 

that is, unless the arguments make it clear, no argument means real space, no argument with a tilde denotes 
Fourier space. 

We can write the Lagrangian halo overdensity as an expansion over the large-scale (p modes, 



Sh= d k Mk) + ^ dk,d k2^^ 



(pe{ki) (j)e{k2) 



dmJ_lad5s'no{Ss,(T'^,Sc) dmJ_''^dSsIlo{6s,a'^,Sc) 

(69) 



where (. . .)o means evaluating at (f)e = 0. Note that we expand only in the large-scale Gaussian field (pe{k). 

We shall sec soon that the perturbation expansion reorganizes itself into local terms in Sf (x) (which contains 
PNG) and non-local terms in (j)i{x). The former are the same terms present for Gaussian initial conditions 
that lend itself to a local (Eulerian or Lagrangian) bias description, while the latter become local in (/)f(x) 
only for local PNG. Therefore, we do not assvme the nature of bias (local or nonlocal, and in which fields), 
it all follows from applying the PBS to the particular PNG model under consideration, which in our case 
they are all described in terms of a Gaussian field ^ and its interactions (PNG kernels). 
The first derivative required in Eq. (69) is, to all orders in PNG 

\V(l)e{k) ) ~ ^^\dJ^ ) oKVcPeik) ) ~ ^Mg) o\V(l)e{k) ) o'^ Qc^p} VD^^(fe)jo ^ ' 



where to first order in /nl, as discussed above, only the first two terms contribute {p = 1,2) while p = 3 
contributes to ©(/nl) ^^'^ C(ffNL), see Eqs. (56) and (61). The p = 1 contribution is the usual PBS bias 
present in the Gaussian initial conditions case, while the rest are new contributions due to PNG. Similarly, 
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for the second derivative we have, to all orders in PNG 



(71) 



VS^2 Jo \V(l)e{ki) ) o\ 



vSi \ /an 



/ a^n \ r/ \ / i^o" \ / Udg \ / Ua^ \ 



^dSedal^Jo l\'D(l)i{ki)J o\'D4>e{k2) 



Vcl)e{k2))o\VMki)^o 



^ / N (^^\ + SIIo / W X 

^ a(a^)2 V©.^^(fei);oV©.^^(fe2)/o^aa^ V©<A^(fei)©.^^(A;2)>'o^"'" ^ ^ 



where in the second equality we singled out the most 
dominant contributions: again the first term here 
is present even for Gaussian initial conditions, the 
second and third are ©(/nl)) and the last line is 
J and 0(5nl). 
We now have all the ingredients to write down 
general expressions for the PBS predicted linear and 
quadratic bias to all orders in PNG. 



D. Linear Bias 

Let us first look at the linear bias terms, and in- 
clude for now only the first derivative contributions 
in Eq. (69). Wc are interested in the bias in Fourier 
space, which is easy to obtain from Eq. (69) since 
the only a;-dependence comes from the first deriva- 
tives of n which from Eqs. (70) and (66) means that 
they are just proportional to a plane wave. Thus we 
obtain for the Fourier space Lagrangian halo pertur- 
bation, 



~6l: 



/pi/(an/ac(f))o 



dm /no 



(73) 



where to simplify notation we omit arguments which 
are displayed in Eq. (69). The sum over p here in- 
cludes all orders in PNG through /pi, the first term 
is the usual scale- independent linear Lagrangian bias 
obtained from derivatives of the conditional mass 



function 

p=l: 



dmj{dIi/d5i)o 



dm /no 



dSf 



In 



f dn{5e) 
\dh\m 



(74) 



whereas the next gives the possibly scale-dependent 
bias correction to C(/nl), 



p = 2: 



,(2) _ a„[j2i/ano/ac7^ 



M(fc)a„/ no 

dal^ [hi ^o] 



(75) 



where in this last equality it is understood that any 
m-dependence inside the square brackets is rewritten 
in terms of cr^. Going from the first-upcrossing rate 
J-Q to the mass function we can rewrite this as, 



1,(2) 



d„ 



'21 



(76) 



This is our final prediction for the possibly scale- 
dependent bias correction due to PNG at leading or- 
der in /nl in terms of the mass function (dn/dlnm) 
and the variance smoothed at the mass-scale of the 
halo CT^ (see Eq. 50 for definition of M). Note that 
to leading order, since /21 depends on the primordial 
bispectrum (Eq. 66) 



l2i{k,m) = 



B^^iq,k-q,-k)d\ (77) 
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the mass function can be computed from Gaussian 
initial conditions, therefore Eq. (76) can be com- 
puted given a PNG model from Gaussian simula- 
tions alone /or any type o/PNG, without making any 
assumptions about universality of the mass function. 
Markovian behavior or having a specific implemen- 
tation of the excursion-set approach. Note that our 
result, Eq. (76), is different from that in [9] where 
the scale-dependent bias for local PNG without as- 
suming universality is given in terms of derivatives of 
the mass function with rosp(X't to (Ts, as that woiild 
require running Gaussian simulations for different 
normalizations. Our expression is for valid for gen- 
eral PNG (characterized by Eq. 77) rather than local 
and one only needs Gaussian simulations with the 
same erg to make predictions. 

That Eq. (75) can be written in terms of the mass 
function without any assumptions about Markovian 
evolution and universality is only because this cor- 
rection comes from a modulation of the variance 
and thus involves dU.Q / da^ which gives the first up- 
crossing rate, and thus can be written back in terms 
of the mass function. This fact is a special property 
of /nl-PNG, on the other hand, at cubic order 



p = 3 



6(3) = 



(3)i 



M{k)dmJUo 



(78) 



and there is no simple way of proceeding further 
without knowing how flo depends on the third mo- 
ment of small-scale primordial fluctuations. Because 
/si depends on m, Eq. (78) cannot be written solely 
in terms of derivatives of the mass function with 
respect to the third cumulant. This is important 
for S'NL-PNG, for which this term makes the leading 
contribution to possibly scale-dependent bias. And 
note that Gaussian simulations are not useful in this 
case, one must resort to some kind of approximation, 
e.g. assuming that IIo can be replaced by the PDF 
of the density field (violating the boundary condi- 
tion that Ila{Sc) = 0) and then (expanding the PDF 
in Edgeworth series to compute its derivative with 
respect to c^^' (see [62-64] for work along these lines) 
or use results from the excursion-set approach in [65] 
to compute how Ho depends on . We won't com- 
ment on this further, except to note that clearly in 
this case computing the scale-dependent bias contri- 
bution relies on a lot more assumptions than for /nl- 
PNG where a clean result such as Eq. (76) is possible 



for arbitrary primordial bispectra. For a calculation 
of scalci-dcpcmdent bias up to trispectrum contribu- 
tions from local bias models see [66]. 

So far to calculate linear bias we have relied on 
Eq. (73) which only includes first derivative contri- 
butions from Eq. (69). In going from Eq. (73) to 
Eq. (74) we have implicitly assumed that Iii4>e = 
M(f)£ = Se which is only true in the Gaussian case. 
However, it's easy to sec that the second deriva- 
tive contributions with p = 1 in Eq. (72) contribute 
the non-Gaussian part of with the same bias co- 
efficient, and thus Eq. (74) holds when comparing 
the halo to the fully non-Gaussian matter perturba- 
tions (higher-order PNG terms inside Si get gener- 
ated by higher-derivatives). This goes to show that, 
as stated in Section IIIC, the perturbation expan- 
sion can be reorganized in terms of (j)£ and Si. The 
remaining contributions coming from second deriva- 
tives contribute to quadratic bias, and will be con- 
sidered below in Section IIIG. 

Summarizing, we can write the Eulerian linear 
bias in Fourier space as 

bi = l + biL 



(2) 



where b^^^ = 1 



,(P) 



(79) 



b'fl (p > 2), and the 



first two terms b^^^ ^^d 6^^' can be computed from 
a given conditional and unconditional mass function 
(Eqs. 74 and 76), but the additional contributions 
from cubic and higher-order PNG cannot be simply 
related to mass functions without further assump- 
tions. 

From now on, we will concentrate on /nl PNG to 
leading order in /nlj and thus the contribution to 
possibly scale-dependent bias (in either Lagrangian 
or Eulerian space) is 



Abi = Mil = b' 



(2) ^ ;,(3) 



IL 



l(2) 



(80) 



thus henceforth we will use A5i to denote the scale- 
dependent bias from /nl PNG. Table I summarizes 
our notation for bias parameters. 
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TABLE I. Different bias parameters used in the text 



Symbol 


Meaning 


First appearance 


L 

01 


linear Eulerian bias (includes scale-independent and dependent contributions) 


Hiq. (/9j 


biL 


linear Lagrangian bias (includes scale-mdependent and dependent contributions) 


Eq. (79) 


Oil 


scale-independent contribution to linear Lagrangian bias 


Eq. (74) 


t{2) 

Kl 


scale-dependent correction to biL from variance modulation ©(/nl) 


Eq. (75) 


.(3) 


scale-dependent correction to biL from third-moment modulation 0{gisi^) 


Eq. (78) 


Abi 


sum of all contributions to scale-dependent bi {Abi ~ b^^l for /nl PNG) 


Eq. (80) 


AbiL 


sum of all contributions to scalc-dcpcndent biL, Abit = Afei 


Eq. (80) 


ba 


bias measured in Gaussian simulations from halo-matter spectrum 


Fig. 5 


bres 


residual halo bias in simulations after theoretical scale-dependent bias is substracted 


Eq. (Ill) 


SblL 


PNG correction to the scale-independent linear Lagrangian bias 


Eq. (113) 


b2 


quadratic Eulerian bias (includes scale-independent and dependent contributions) 


Eq. (116) 


b2L 


quadratic Lagrangian bias (includes scale-independent and dependent contributions) 


Eq. (104) 




scale-independent contribution to quadratic Lagrangian bias 


Eq. (93) 


=,(1,2) 


scale-dependent correction to 62L from mean and variance modulation ©(/nl) 


Eq. (100) 


J,(2.2) 
"2L 


scale-dependent correction to b2L from variance-variance modulation ©(/nl) 


Eq. (101) 


"2L 


scale-dependent correction to 62L from variance modulation ©(^nl) 


Eq. (102) 



E. Comparison with Known Results 



function 



Here we compare our results to known PBS re- 
sults in the literature, for comparison with other 
than PBS approaches see Section V. 

Let us first evaluate these results for local PNG of 
/nl type, i.e. K = 1 and no cubic or higher-order 
PNG. To make contact with the literature, we work 
in the low-k limit where (sec Eqs. 63 and 77) 

7^i^(fc ^ 0, m) = 4/nl (tI + 0{k^) (81) 



and then Eq. (75) reads 



A6i(fc ^ 0) 



4/1 



NL 



M{k) 



ai„,2 In(a^J-o) (82) 



where again as in Eq. (75) any dependence of J-q 
on m is rewritten in terms of a^. This result is 
still more general than those in the literature, which 
assume Markovian behavior and universality of the 
mass function to relate the scale-dependent bias am- 
plitude to the scale-independent Lagrangian bias. 
To obtain this limit, let's assume Markovianity, in 
which case the PBS linear Lagrangian bias can be 
written as a derivative of the unconditional mass 



6(1) - _^ 



and in addition assuming universality (see Eq. 40) 
we can rewrite this as 



(84) 



which reduces Eq. (82) to the well-known result [9, 
10] 



M{k) 



(85) 



In section IV we will compare N-Body simulations 
to this standard result and our more general result, 
Eq. (76). The latter incorporates corrections to this 
formula due to three effects: i) non-universality of 
the mass function, ii) non-Markovian behavior, iii) 
beyond-leading order corrections in k to Eq. (81) 
coming from M (see Eq. 55), the transfer function 
and the smoothing kernel. These lead to correc- 
tions to the scale-independent bias to ©(/nl), see 
Appendix D for more details. 
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FIG. 2. Full calculation of hi, Eq. (77), for local PNG 
normalized by its large-scale limit, Eq. (81), for halos of 
mass logj^o = 16,15,14,13,12 (dashed, dotted, solid, 
long-dashed, dot-dashed). 



FIG. 3. Same as Fig. 2 but for non-local PNG of orthogo- 
nal (top panel) and equilateral type (bottom panel) . The 
J21 are normalized by the large-scale values predicted by 
Eqs. (87) and (89), respectively. 



To appreciate this latter point. In Figure 2 we 
show the full calculation of /2i(fc, m) as a function of 
k for different halo masses, normalized by it's large- 
scale limit, Eq. (81). We use this numerical calcu- 
lation of I21 in all our computations below, except 
in Appendix D where the different sources that con- 
tribute to scale-independent bias are separated by 
Taylor expansion. 

For non-local PNG, assuming Markovianity and 
universality we can rewrite Eq. (76) or its equivalent 
Eq. (75) as 



A61 



hi{k,m) (1) 1 
2alMk) ' M{k) 



9ln era. 



hi{k,m) 



(86) 

The first term here agrees with the results of [33], 
the second term (which vanishes for local models 
in the low-fc limit) has been recently taken into ac- 
count in [62, 63]. Again, our results in Eq. (76) are 
more general than these as they incorporate the ef- 
fects of non-universality of the mass function and 
non-Markovian behavior. Eq. (86) is of course nice 
in that it's analytic (aside from the computation of 
/21), whereas Eq. (76) requires knowing the mass 



function. However, let us re-emphasize Eq. (76) re- 
quires only the Gaussian mass function, and thus 
accurate fitting formula could conceivably be used, 
though in this paper we use the measured Gaussian 
mass function (see Eq. 112 below for details of how 
we implement this). 

Let us now specifically evaluate these results for 
the orthogonal and equilateral templates, in partic- 
ular their scale dependence, which is governed by 
/21. In the low-A: limit, we can expand the cross 
bispectrum in powers of {kg /kg) with ki ^ fc^, the 
so-called squeezed limit, which is described in de- 
tail in Appendix E. For the orthogonal template, 
using Eq. (E16) for the squeezed limit of the cross 
bispectrum we obtain, 

I^{\k ^ 0,m) « -12/NLfcSi(m), (87) 

where we introduced (note that So(m) = cr^) 



S„(to)= / d\P{q)W^^{qR) q- 



Equation (87), when used into Eq. (86), gives the 
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expected k^^ scale-dependent bias^, while for the 
equilateral model, using Eq. (Ell) for the squeezed 
limit of the cross bispectrum we obtain instead, 



(89) 



which leads to fc"*"^ scale-dependent bias, although 
with a very small amplitude (percent level for /nl = 
—400, see Fig. 8 below). This is not the expected 
result, but it arises because of cancellations when 
integrating over the cross bispectrum, while a naive 
use of the scaling of the kernel would predict a 
scale-dependent bias, see Appendix E for details of 
the calculation. This different scaling is made possi- 
ble in the first place because we used (/> instead of $ 
as the split variable, and for the equilateral template 
short ^-modes can contribute to $ at large distances 
through the singular squeezed limit of the K^^^ ker- 
nel. See the next subsection and Appendix E for 
more discussion. We will compare both approaches 
with simulations in Fig. 8 below. 

Figure 3 shows the results of the exact numerical 
integration of I21 for these non-local models com- 
pared to their large-scale limits. We see that there 
are next-lo-leading order corrections to the large- 



scale result similar to the local PNG case shown in 
Fig. 2, that are particularly strong for the equilateral 
case (bottom panel, where drops to zero faster 
than in the other PNG models). 



F. Using 3> as the large-scale constraint 

Let us now go back and ask what would have hap- 
pened if we imposed the large-scale constraint on $ 

rather than (p. The procedure is the same, but the 
expansion of the halo perturbation in Eq. (69) is now 
in ^i(k) rather than 4>e{k) and we need to calculate 
the modulation of the cumulants of the small-scale 
density field 6s by the large-scale field which 
is non-Gaussian. This modulation is described by 
the general result in Eq. (66), of which Eq. (64) is 
the simplest version, corresponding to modulation 
of the variance by the large-scale field at 0(/nl)- 
These results are exact for 0^ a Gaussian field, as 
one can easily check by doing a Taylor expansion of 
c^P"> on and computing left- and right-hand sides. 
We need to extend these results to the non-Gaussian 
case when (p^ is replaced by 

Let us consider the simplest case first, Eq. (64). It 
is easy to see that in the non-Gaussian case it gets 
modified to, 



(—) 



P^{k)P^{q)P^{\k-q\) 



B^{q, k- q, -k) + ... 



(90) 



where we neglected higher-order terms in PNG. The first term leads, to leading order, to the same Eq. (64) 
but where -Bjj^ is replaced by B-g^^, while the second term is at least of 0(/nl) can be dropped. This 
implies that the scaling of the halo bias in the case of the equilateral template will be changed from k~^^ to 
scale-independent. See Fig. 8 below for comparison of this prediction against simulations. For the other two 
templates (local and orthogonal) there are no significant changes in the predictions of the scale-dependent 
linear bias. 

The situation is a bit more interesting for the case of quadratic bias that will be discussed in the next 
subsection. For this we need in addition the modulation of the variance by two large-scale fields, which 



^ Recall, however, as mentioned in the introduction, that this 
is due to a pathology of the template we use, the true or- 



thogonal PNG should scale as 721 
independent bias instead. 



k leading to scale- 
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reads^ 

( T^^CJ-" N ^ (a^<&Kfcl)'&|(fc2))c B^{ki2,-kl,-k2) ((7^$;(fci2))c . . 

VP$^(fcl)P$^(fe2)yo P^{h)P^{k2) P$(fcl)P*(fc2) P$(fcl2) ^ ^ 

where now the first two terms are of the same order, i.e. C(/nl) for /nl-PNG. These results can be cast in 
terms of the /pg's defined in Eq. (66), by saying that when modes have been constrained we must change, 

/2l(fc) ^ hl{k), /22(fel,fe2) ^ /22(fel,fe2) - /2l(fcl2) X ^"i^'f ' 7^'"^'' (92) 

P$(A:i)P$(fc2) 

where the Ip^'s are calculated in the same way as Eq. (66) but with T^t'j ^ replaced by r>t«^ ^. We 
now discuss the implications of these results for quadratic bias. 

G. Quadratic Bias 

Let us now go back to Eqs. (69-72) and collect the terms which lead to quadratic Lagrangian bias. Of the 
five terms displayed in Eq. (72), one of them (the second on the first line) was already taken into account 
as it contributes to linear bias (restoring the non-Gaussian part of Se). The remaining four contribute to 
quadratic bias. Because of Eq. (66) all terms come with an a;-dependence that is just a plane wave with 
momentum ki2 = ki + k2, and thus the Fourier space biases are straightforward to compute. The first term, 
also present for Gaussian initial conditions, leads to the scale-independent quadratic bias of local form 

that is, it contributes a term b^^6j/2 in the expansion of 6h given by Eq. (69), or in Fourier space 

biti^l^e, 5e] = b^^l j [6j,] 5i{k^) Seik^) d^k,d^k2, (94) 

where [Sb] = S^ik — fei2). Recall that Si is the non-Gaussian large-scale density perturbation. While this 
contribution from Eq. (72) leads to only the Gaussian part of Sg, third and fourth derivative terms in Eq. (69) 
with the same coefficients restore the non-Gaussian parts, as it is easy to check. 

The other three terms are not local in real space (except for local PNG, see below), and they generically 
depend on scale. From the second line in Eq. (72) we obtain to ©(/nl) 



d„ 



rUm [l2i{ki,m)+l2i{k2,m)]J{d^n/d6eda'^)o 
2LkiH^^^Si]= / — a f TT ^-[6Ti](j)i{ki)Se{k2)d^kid^k2 (95) 

J elm j i-lo 



while the other contributions are C(/nl) 



drr 



/2i(fci, m)/2i(fc2, m) / d^no/d{alf 



5m /no 



[5T,](t>t{ki)(^e{k2)d''kid'k2 (96) 



5 These ealeulations are precisely equivalent to those in RPT ^^^^^^ ^j^^;^ ^q. (33), while Eq. (90) fixes a sign typo in 

that relate multi-point propagators to cross-correlations, their Eq (34) 

see [67] for the non-Gaussian case. Equation (91) corre- 
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and O(ffNL) (and also 0(/nl)) 



J2) 



[<i>e,'Pe] 



d„ 



I22{ki,k2,m) JdUo/da^ 



dm 1^0 

Putting all these together in Fourier space we have 



>1L 



(97) 



(98) 



The superscripts in 62L denote the physical origin of 
each term, e.g. (1,2) in Eq. (95) indicates that this 
contribution comes from modulation of the large- 
scale mean and small-scale variance, sec also Table I 
for explanation of our bias parameter notation. Note 
that while is a Gaussian field, includes PNG 
contributions (which are higher-order in PNG be- 
cause of the overall factors of /pg). If we ignore these 
higher-order contributions, we can simplify the no- 
tation by defining quadratic bias kernels in Fourier 
space in terms of the J^'s, e.g. 



,^t\ = J bif\ki,k2)ST,{k-ki2) 

X 6i{ki)Si{k2)(fki(fk2, (99) 



and similarly for the other two cases. The leading- 
order contribution 0(/nl) to scale-dependence of 
quadratic bias is thus given by Eq. (95), which from 
Eq. (74) can be rewritten in a compact form (com- 
pare to Eq. 75) 



biY\kuk2) 



d,2jl2liki)b['lj'0 

M{ki)To 



+ ki^k2 



(100) 

where again it is understood that any mass depen- 
dence inside the square brackets is rewritten in terms 
of fT^, and we have suppressed the mass arguments 

of /21, s-nd for simplicity. The other two 
kernels can be similarly written, 



,(2,2) n 1 ^ S< [^21 (^l) ^21 (fc2) 9<^^^o] 



^^(fci.fca) 



M{ki)M{k2)J'o 



8^2^ [I22{ki,k2) To] 
M(fci)M(fc2)J-o 



(101) 



(102) 



These results can also be cast in terms of the 
mass function, in analogy to the linear bias scale- 
dependence, Eq. (76). For example, for the leading 
correction to C(/nl) we have from Eq. (100) 



b^^i^\kuk2) = 



d„ 



l2likl)b['l 



dn 
d\nn 



y dm J 



M{ki) (^) 
ki o k2 



(103) 



which again can be implemented using measure- 
ments in Gaussian simulations as discussed above. 
Summarizing, the total quadratic bias kernel is sim- 
ply the sum over all these contributions, i.e. 



'J2L 



^(l) ^ ^(1.2) 



"2L 



^2L 



f,(2,2) 
^2L 



"2L 



(104) 



where only the first term is of local form for non- 
local PNG. 

Following the treatment in Section HIE, wo can 
take the low-A; limit and assume universality and 
Markovian evolution to simplify these expressions, 
and make contact with previous literature. For local 
PNG we have using Eq. (81) that the leading-order 
scale-dependence, Eq. (100), reads 



(,(1.2) _ r, . /c ,(1) 

O2L ~^/nl(c>c02L 



1 



1 



M(fci) M{k2) 



while for Eq. (101) we have 



^(2,2) 
"21, 



4/: 



M{ki)M{k2) 



(105) 



(106) 



To give the low-fc limit of Eq. (102) we need (from 
Eqs. 56-58 and Eq. 66) 

I22{kuk2 ^ 0, m) = (8/^L + 125nl) <7™ (107) 
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which gives (see Eq. 46 for definition of 5nl) 



2,(2) 



(1) 
IL 



M(fci)M(fc2) 



(108) 



Equations (105-106) and (108) agree with previous 
results in the literature for local PNG [29, 68]. There 
are no results on non-local PNG quadratic bias in 
the literature. 

Let us now discuss how these results change if 
we impose the large-scale constraint on rather 
than 4)1. It's easy to check that the same results 
in Eqs. (100-103) hold for quadratic biases, but with 
the /pq's changed by /pq's according to Eq. (92). For 
local PNG, is still true that Eq. (81) holds in the 
new variables, i.e. = 4/nl cr.m, but Eq. (107) 
gets changed to 



+ ki ^ k2 

+ (8/SL + 6gNL) 



P<t>^\ki-q\) 



d\P^{q)P^{q) + 12 g^l^ a, 



2 



(109) 



where Pj(<z) is the power spectrum of the smoothed 
density field. In the low-A; limit, the factor in 
square brackets vanishes, therefore, only the first two 
terms give a modification of the ^^(/nl) amplitude 
in Eq. (107), while the (7nl amplitude remains the 
same in this limit. Therefore, we see that Eqs. (105- 
106) are unchanged but Eq. (108) now reads. 
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This modification of the ©(/nl) quadratic bias pa- 
rameter may be probed through measurements of 
the halo bispectrum as a function of triangle shape. 
It arises from the same effect that can change the 
scaling in the low-fc limit for the linear bias, from 
contributions of the K'-^^ kernel that couples two 
short 4>s modes. 

Finally, note that these results are for the La- 
grangian quadratic bias parameters, what we need 
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FIG. 4. The bias for FOF0.156 halos as a function of 
scale for Gaussian and local, orthogonal and equilateral 
PNG initial conditions. Since cosmic variance is domi- 
nated by its Gaussian contribution, we only show error 
bars on the local PNG case for clarity. The equilateral 
and Gaussian case are very close to each other, whereas 
the orthogonal template (/nl = —400, blue dashed lines) 
is in between them and local (/nl ~ 100, red dotted) for 
significantly biased objects (top two panels), but below 
the Gaussian (black solid) and equilateral (/nl = —400, 
green solid) case for low-mass halos at 2; = 0. 



to compare against simulations is to compute their 
Eulcrian counterparts. This is a standard procedure 
usually done in the spherical collapse approximation 
(see e.g. [29, 30, 68]) or, more accurately, full per- 
turbation theory. We leave this for an upcoming 
work where we implement these PBS predictions for 
the bispectrum and compare against simulations for 
halos and mock galaxy catalogs. 



IV. COMPARISON WITH SIMULATIONS 

We now contrast our predictions for large-scale 
linear bias with measurements in the simulations 
discussed in section HE. Since our predictions for 
the scale-dependent bias from PNG should be more 
widely valid than the standard results based on uni- 
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versality and Markovian evolution, our primary goal 
here is to test for the amplitude of this scale de- 
pendence. Previous results in the literature on this 
proceed by modeling the full bias factor, including 
scale-dependent and independent contributions, see 
e.g. [40, 62, 63, 68-72], and there is no consensus 
about whether a 'fudge factor' is needed to properly 
account for the amplitude of scale-dependence for 
local PNG. 

There are many reasons why this might be the 

case. First, not all works used the same halo def- 
initions, we explore the dependence on halo defi- 
nition below. Second, there is the impact on halo 
bias from from transients induced by setting up ini- 
tial conditions in the simulations [41]. For exam- 
ple, we find that using Zel'dovich initial conditions 
instead of 2LPT at z = 49 for local PNG with 
/nl = 100 leads to a 2: = 1 halo power spectrum 
{M = 1013 - IQi^'Mo/Zi) that is larger by 14% at 
k = 0.003 /iMpc"i and 3% at A: > 0.05/iMpc"\ 
These transients also induce artificial violations of 
universality. 

From the theoretical point of view, deviations 
from the standard predictions are expected by vi- 
olations of Markovianity and universality. While 
deviations from the former have not yet been es- 
tablished in a precise quantitative way, there is a 
significant body of work showing that universality 
of the mass function does not hold at the 5-10% 
level [41, 53, 54, 73-75] for FOF halos, with more sig- 
nificant deviations for spherical overdensity (SO) ha- 
los [75]. In addition, the peak-background split cal- 
culations for Gaussian initial conditions show sim- 
ilar deviations [53, 54, 76]. In this case, however, 
there is the extra complication in going from the 
bias parameters in the expansion of perturbations 
to the bias parameters that appear in the correla- 
tors such as the power spectrum, which will differ in 
general by renormalizations induced by loop correc- 
tions [52, 77, 78]. 

Figure 4 shows the bias computed from the halo- 
matter power spectrum for one of our choices of halo 
definition (FOF halos with linking length equal to 
0.156 times the interparticlc separation) as a func- 
tion of scale for Gaussian and local, orthogonal and 
equilateral PNG initial conditions. We see the ex- 
pected scale dependence for the local case, a weaker 
dependence for the orthogonal template, and close to 
Gaussian bias in the equilateral model. Rather than 



performing a global fit for the scale-independent and 
dependent terms, our approach here is to look at 
the residual halo bias in simulations after the scale- 
dependent bias predicted by theory is substracted, 
i.e. (see Eq. 80 for definition of A61) 



^ -'mass ' Nbody 
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where Phm is the cross-spectrum between halos and 

matter. Note that the N-body quantities arc for the 
PNG model imdcr consideration, i.e. the mass power 
spectrum includes PNG. Simulations and perturba- 
tion theory calculations show that there are interest- 
ing PNG corrections for the mass power spectrum 
and bispectrum (see [67, 79, 80] and Fig. 1), and 
even down to the nonlinear regime [81], but we won't 
explore those here. 

If the theoretical model Abf^°°'^^ is correct, the 
residual bias 6res should be consistent with scale in- 
dependence, whereas if the theoretical model does 
not predict the correct scale-dependent bias feres 
will still show residual scale-dependence. Further- 
more, provided that feres is consistent with scale- 
independence, we can look at the ratio of feres to be, 
the halo bias measured in our Gaussian simulations, 
to quantify the magnitude of the PNG corrections 
to the scale-independent bias. 

To calculate our predictions for Afe'^'"^°'^'^, given by 
Eq. (76), a mass derivative of the Gaussian mass 
function is required. We implement this by doing. 
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where Nh{i) denote the number of halos in the 
Gaussian realizations in a bin of constant din to, 
Nf°^ = J2i Nh{i), and the sum is over the mass bins 
belonging to the halo sample. The numerical deriva- 
tive (denoted by a prime) is taken by doing cen- 
tral differences from neighboring bins. Care must be 
taken at low mass to have a smooth mass function, 
when the number of particles in a halo is smaller and 
binning effects can indiice artificial noise. 

In Figure 5, we show feres for FOF halos with 
linking length equal to 0.2 times the mean inter- 
particle separation with local PNG with /nl = 100 
normalized by the halo bias measured in our Gaus- 
sian simulations 6g (as labeled in each panel), for 
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FIG. 5. The residual halo bias in local PNG with 
/nl = 100 (normalized by the bias measured in Gaus- 
sian simulations be) after the theoretical PBS scale- 
dependent bias is accounted for using two predictions: 
our result (Eq. 76, blue squares), and the standard pre- 
diction (Eq. 85, red triangles). This is for FOF0.2 ha- 
los and different halo masses and redshifts. Our pre- 
dictions are consistent with scale-independent residuals, 
while the standard prediction is not, more so for large- 
bias objects. The dashed lines show the expected ratio 
of scale-independent biases assuming universality plus 
Markovianity. 



different halo masses and redshifts. The symbols 
with error bars show 6rcs/^G for two different the- 
oretical models, our prediction Eq. (76) shown by 
blue squares, and the standard prediction Eq. (85) 
denoted by red triangles. We see that our predic- 
tion for the scale-dependent bias performs better, as 
the residuals are consistent with scale-independence, 
whereas the standard prediction is not. The latter 
over-predicts the amplitude of the scale-dependent 
bias, as a result the residual bias 6ios is suppressed at 
low-fc. Figure 6 shows the analogous results for FOF 
halos obtained from a linking length 0.156 times 
the mean interparticle separation, and shows a sim- 
ilar overprediction of the scale dependence by the 
standard formula. The magnitude of this deviation 
is somewhat larger for FOF0.2 halos, thus the de- 



FIG. 6. Same as Fig. 5 but for FOF0.156 halos. Similar 
results hold despite a very different halo definition. 



tails depend on halo definition. This is in qualita- 
tive agreement with previous studies that required 
a "fudge factor" less than unity (typically q ~ 0.75) 
on top of the standard prediction [28, 69, 70, 82]. 

From the constancy of the ratio bios/^G for 
our theoretical prediction (square symbols) we can 
read off that there is a PNG correction to scale- 
independent bias. The sign of the magnitude is ex- 
pected as for local PNG with positive /nl the halo 
mass function is enhanced and the scale- independent 
bias is thus suppressed compared to the Gaussian 
/nl = case. To be more specific, we show using 
dashed lines in Figs. 5 and 6 the expected scale- 
independent correction to halo bias assuming uni- 
versality plus Markovianity, that is (see Eq. 84) 



db' 



(1) 

IL 



where 



u(l)| 



^^^^ ((5r)( dm ) ™ (dlnm) 



^ilIg 



(113) 



^dlnm/ PNG 

(114) 

and similarly for the Gaussian case [9, 25, 68-70]. 
Because this is for fixed mass, we integrate each ex- 
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pression for the bias at fixed m weighted by the cor- 
responding mass function over the desired mass bin. 
We see from Figs. 5 and 6 that these predictions, 
for a wide set of halos (note the range in Gaussian 
bias parameters from 1.38 to 6.23) match rather well 
the residual bias from our theoretical prediction, al- 
though there are certainly deviations at the percent 
level. This fact, together with the flatness of the 
residual bias as a function scale, tells us that our 
improved treatment leads to a better description of 
the amplitude of scale-dependence in local PNG. 

In [71] it is found that spherical overdensity (SO) 
halos obey the standard formula for scale-dependent 
bias more closely than FOF halos. Naively, since 
SO halos violate universality more strongly than 
FOF halos [75], one would have expected the oppo- 
site (particularly at low mass, where the deviations 
from PBS bias plus Markovian and universality are 
stronger [76]). We don't currently have SO halos for 
the simulations we present here, but would be inter- 
esting to check our improved theoretical prediction 
against SO halos. 

Figure 7 shows the analogous results for the or- 
thogonal template with /nl = —400. Our predic- 
tion for the amplitude of the bias (blue squares) 
leads again to a residual bias consistent with scale- 
independence at low-fc. Note from Fig. 4 that in 
our low-mass bin at z = (bottom panel) the 
scale-dependent bias changes sign, and our predic- 
tions correctly match this (second panel from top in 
Fig. 7). At high-fc, as nonlinear scales begin to be 
probed around k ~ 0.1/iMpc~^ there is significant 
evidence for scale-dependent non-Gaussian contribu- 
tions unlike the local PNG case shown in Fig. 6. This 
must be due to the larger value of /nl in the orthog- 
onal case (/nl = —400 versus /nl = 100 in the local 
case). We also show two other predictions, assuming 
Markovianity and universality (red triangles) which 
gives rise to Eq. (86) (as in [62, 63]) and in green 
pentagons its first term only (corresponding to the 
predictions in [33]). While this is for FOF0.156 ha- 
los, we find very similar results for FOF0.2 halos. 
We conclude that our improved formula performs 
best compared to the alternatives. Note in this case 
that the residual bias predicted by Eq. (114) (shown 
as dashed lines in Fig. 7) shows larger deviations 
than for the local case. This might be due to non- 
Markovian corrections proportional to /nl [61] that 
are not included in Eq. (114). 
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FIG. 7. Residual halo bias for orthogonal PNG with 
/nl = —400 after our predictions for scale-dependent 
bias are included, Eq. (76), for FOF0.156 halos. The 
blue square symbols show our predictions, red triangles 
the predictions of Eq. (86), and green pentagons its first 
term only. For our predictions, residuals are consistent 
with scale-independence at large scales for different halo 
masses, redshift and halo definitions. 



Finally, in Fig. 8 we present residual halo bias re- 
sults for equilateral PNG for FOF0.156 halos. In 
this case we compare our prediction for residual 
bias (blue squares) based on subtracting the scale- 
dependent term given by the bottom panel in Fig. 3 
(and Eq. 89 in the low-fc limit) and without substrac- 
tion (red triangles) which correspond to the standard 
prediction (that includes only scale-independent cor- 
rections) and also the PBS prediction when the con- 
straint is done on the field (see Section lllF). 

We see that the differences are small, although 
the measurements are slightly more consistent with 
scale-independent residual bias for the prediction 
based on the cross-bispectrum -B^*^ rather than 5$, 
but we don't consider this statistically significant. 
Note that the sign of the scale-dependent effect in 
this case depends on halo mass (negative for high 
mass and positive for our low mass bin), and for 
high-mass at z = 0.97 (with Gaussian bias ba — 6) 
the effect is only about 2%, thus for all practical 
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FIG. 8. Residual halo bias for equilateral PNG with 
/nl = —400 after predictions for scale-dependent bias 
are included (blue squares), versus assuming no scale- 
dependent contributions as in the standard predictions 
(red triangles). The different panels show two halo def- 
initions, redshift and halo masses. Note the enhanced 
vertical scale in this figure. 



purposes not very important. 



ON LOCAL BIAS VS PBS, LOOPS AND 
EFFECTIVE THEORY OF BIAS 



Let us now concentrate in this section in the sim- 
plest case, local PNG, where most results in the lit- 
erature have been obtained. We will discuss first the 
large-scale behavior (scale-dependent bias at low-fc) 
which arises from loop corrections in local bias, and 
how this compares with simulations. We then dis- 
cuss what happens with bias towards the non-linear 
regime in the context of loops in the PBS, which 
provides a nice example of a formulation of an effec- 
tive theory of bias that can be applied to the more 
complicated case of galaxies. 

In Eulerian local bias model [83] the smoothed halo 
perturbation Sii{x) is assumed to be a local function 



of the smoothed dark matter perturbation 
Snix) ^ mx)] 



(115) 



Note that both fields are smoothed on the same 
scale, say Rs, which is the scale one is interested 
in measuring things at, and must be large enough 
compared to the Lagrangian size of halos for local 
bias to make sense at all (operationally the function 
/ can be obtained from simulations by constructing 
a scatter plot such as the one presented in Fig. 9 
below). At large smoothing scales, density pertur- 
bations are small so one can write a perturbative 
version of this 
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where we have suppressed a constant term 5o that 
enforces {5h) — 0. This leads to a halo-matter power 
spectrum. 
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where B denotes the bispectrum of the density field, 
and we have deliberately left the UV cutoff in the 
loop momentum unspecified. Interestingly, as noted 
in [24] the primordial contribution to the density 
bispectrum gives rise to a fc~^ scale-dependent bias 
proportional to /nl, i-e. 
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where is the variance of the smoothed density 
field. In such an approach, the sign of Abi is pre- 
cisely given by that of &2- This result seems quali- 
tatively similar to Eq. (85), but the amplitude de- 
pends on 62 and instead. The former, when mea- 
sured from a scatter plot between Sh and S, does not 
depend sensitively on the smoothing scale Rs (for 
Rs > 30 h-^ Mpc [52, 54]), but of course cr^ is very 
sensitive to Rg, and thus the magnitude of this ef- 
fect is uncertain^. However, if we take the rare event 



^ If we had followed the logic of Eq. (116) the smoothing scale 
would be related to the observation scale k, as in Eq. (115). 
See later in this section for more discussion on this. 
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FIG. 9. Left panel: The Sh vs S relationship for local /nl ~ 100 for halos of mass in the 10" - 10" Mq/H range. 
The mean (solid line) is reasonably described by a local relation with bias parameters bi = 1.34 and &2 = —0.58. 
The standard deviation of 5^ given 5 is shown by dashed lines. Right panel: The scale-dependence of the halo bias 
factor. Despite having 62 < 0, these halos show a scale-dependent enhancement of their power spectrum at low-fc, 
rather than the suppression predicted by the local bias model. 



limit, corresponding to halos (or peaks) of very large 
mass, and we take the smoothing length to be the 
Lagrangian radius of the halo Rl, then b2(J^ « b'^l 
and Eq.(118) agrees with Eq. (85) [79]. Note how- 
ever that this requires, effectively, to have taken the 
UV cutoff in the loop in Eq. (117) to be close to 
i?^^ where local bias is a dubious assumption, and 
even if it were to hold including such modes in S 
would likely invalidate the perturbative expansion 
in Eq. (116). 

It is also interesting that for halos away from 
the high-mass limit, where 62 becomes negative, 
Eq.(118) predicts that the halo power spectrum 
should be suppressed rather than enhanced (for 
/nl > 0) at low-fc. In Fig. 9 we test this prediction 
for FOF0.2 halos of mass in the lO^^ - lO" Mq/H 
range at z = 0. The left panel shows a scatter plot 
of the halo versus matter perturbation smoothed on 
Rs = 50/i~^Mpc scales, showing that these halos 
have &2 < due to the negative curvature of the 
solid line, which denotes the mean of the relation. 



The right panel shows however that the bias (de- 
fined from the cross spectrum) is enhanced rather 
than suppressed. Note that another way to see that 
62 < for these halos is that b{k) decreases with 
scale for k > 0.1 /i Mpc~^, as expected since for 
Gaussian initial conditions the sign of 62 controls 
the behavior of the scale-dependent bias at nonlin- 
ear scales [52]. One may argue that 62 calculated at 
Rs ~ 50/i~^Mpc is not fair, since as we discussed 
Eq. (118) approximately agrees with the PBS result 
only when cr^ is taken to be smoothed at the La- 
grangian radius of the halo. However, evaluating 62 
at such small scale does not change the conclusions, 
as 62 is driven to small values (but still negative) 
as the smoothing length decreases [54] . In addition, 
one can see from scatter plots at such smoothing 
scales that local bias is not really a good represen- 
tation of the data. As we discuss below, the PBS 
tells us that we shouldn't push the smoothing scale 
to the Lagrangian scale of halos as the perturbative 
expansion breaks down. 

A generalization of Eq. (118) is available in other 
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flavors of local biasing. In [11], a general formula 
is given for biasing of exponential form, which re- 
produces Eq. (118) to lowest order. Similarly, in 
local Lagrangian bias, [84] derive the same formula 
as [11], and Eq. (118) follows with 62 replaced by its 
Lagrangian version Ij^l for clustering in Lagrangian 
space (see e.g. [33] for recent discussion on this). 
For high-peaks, the Lagrangian to Eulerian map- 
ping does not matter (it is subleading in v'^), but 
for the halos considered in Fig. 9 the contributions 
from the mapping are not negligible. If the mapping 
is done in the usual spherical approximation, then 
the same Eq. (118) follows for the Eulerian (as mea- 
sured in the simulations) scale-dependent bias where 

62 = b^2L + (8/21)^iL) which means that 62L is even 
more negative than the measured 62. In any case, the 
conclusion is that these local bias models (at least in 
their perturbative form) do not match the halo clus- 
tering behavior seen in Fig. 9 where the (Eulerian 
linear) bias is enhanced at low-fc while the quadratic 
Eulerian bias is negative. It would be interesting to 
do a more rigorous peak-type calculation that goes 
beyond a perturbative expansion to see whether the 
predictions differ from a quadratic bias effect for low 
bias tracers. 

Another proposal for dealing with bias in PNG 
models is presented in [85], based on earlier work 
for Gaussian initial conditions [78]. In this case the 
starting point is the Eulerian local bias expansion 
in Eq. (116), but rather than being restricted to 
smoothed fields, 5 is formally assumed to be un- 
smoothed, leading to a scale-dependent bias at low-fc 
for local PNG from Eq. (118) that is formally ultravi- 
olet (UV) divergent. To cure this divergence, a new 
(non-local in density) term proportional to cf) with 
arbitrary amplitude is added to Eq. (116). Since 
this has the same effect in the low-fc power spectrum 
(leading to a fc~^ bias) one can tune its amplitude 
to cancel the UV divergence from Eq. (118). In this 
way, the model is left with a "rcnormalized" finite 
amplitude of scale-dependent bias that can be fit to 
simulations. In a sense, in this approach the UV 
sensitivities arc used as a guide to select the appro- 
priate modification of the local bias model for a given 
PNG type, and after the renormalization procedure 
is complete the theory is no longer UV sensitive. 
By contrast, in the PBS calculation we presented in 
this paper, the effective large-scale degrees of free- 



dom appear naturally after implementing the split 
and there are UV sensitivities left (although prop- 
erly suppressed by inverse powers of the UV cutoff 
of the theory, see below), which alert us to when the 
perturbative expansion breaks down. 

In the PBS, there is also a quadratic bias param- 
eter (see Section IIIG). Does this mean that one 
gets, in addition to the already discussed PBS re- 
sults, another contribution to scale-dependent bias 
from the effect leading to Eq. (118)? The answer 
is 'yes', but since the expansion in the PBS is in 
the large-scale modes dg (with aj <C cr^) the contri- 
butions from quadratic bias are highly suppressed. 
That is, the contribution from the effect leading to 
Eq. (118) reads in PBS (compare to Eq. 85) 
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where for simplicity in the last equality we have as- 
sumed the high-peak limit. Therefore, we see that 
as long as <C cr^ these extra corrections are well 
under control. Since this correction arises as a loop 
integration over the momenta of large-scale modes, 
here is whore the details of the split enter. In our 
formulation, we take the split right at the scale fc we 
are computing the power spectrum (i.e. we choose 
a sharp-fc filter, as explained in Eqs. 47-54), thus 
al is fc-dependent (see Eq. 48). For example, at 
fc = 0.01/iMpc-\ for a 10^^ (10^^) Mq/H halo, 
(c^f/f^m) is 0.0005 (0.0015). At fc = 0.03/7,Mpc-\ 
where the scale-dependent contributions are already 
subleading (see Fig. 4), the large-scale to small-scale 
variance ratios become 0.015 and 0.045, respectively. 
That loops lead to cr|(fc) factors was already used in 
[86] to implement the perturbative bias expansion 
for PNG models of type. 

This shows that the PBS calculation we have pre- 
sented before is self-consistent, i.e. insensitive to 
the inclusion of non-linear bias loops (they vanish 
as fc — >■ 0). In fact, as (ff/cr^) becomes non- 
negligible, one should go back and inchide such cor- 
rections in Eq. (43), which correspond to changing 
cr^ — )• a'^ — af for a sharp-fc filter (see Eq. 56), 
well-known from conditional mass functions [47, 87] . 
In this regime, however, one is trying to study bias 
at scales comparable to the Lagrangian size of the 
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halo, thus finite-size efl:ccts become important and 
bias ceases to be local (in the UV sense, i.e. terms 
with positive powers of k show up). Indeed as k 
approaches the inverse Lagrangian size Rl of halos, 
((T^/(T^) cx (fci?^)"^'^, where n is the effective spec- 
tral index at k. At these scales, the perturbative 
expansion of PBS breaks down, and it is replaced by 
a non-perturbativc, non-local description with ex- 
clusion effects that drive the halo autocorrelation 
function to —1, see [49] (physically, what happens 
is that since Sf is approaching Sc the region of inter- 
est will correspond to a halo, thus proper account of 
the 'cloud in cloud' problem is essential). In the case 
of galaxies, it is at these scales where the non-local 
bias due to halo profiles (in the language of the halo 
model) takes over. 

In [68] a bivariate expansion of the halo bias in 
terms of the smoothed (f) and 5 is developed and it is 
argued that the quadratic effect must be small based 
on choosing a smoothing scale R = 10 Mpc. As 
we mentioned in Section III C this type of approach 
arises from using a real-space top-hat filter to make 
the split between low-A; and high-fc modes, which is 
more appropriate to calculate counts-in-cells statis- 
tics than Fourier space correlators. Our conclusion 
here is similar, but with a key difference: rather than 
choose an arbitrary fixed smoothing scale R, we in- 
stead use a varying splitting scale between small and 
large-scale modes that follows the fc-mode we are in- 
terested in. Since the Lagrangian halo bias param- 
eters can be written as 6„l = S„(i^)/cr^ for some 
dimensionless B„(z^) [51] the perturbative expansion 
of PBS reads (here we ignore terms depending on (f)g 
since we are interested in the UV limit), 

Sf: = Y.^nL5^ = (— )" (120) 

n=l n=l ^™ 

that is, all loop corrections coming from nonlinear 
bias (such as Eq. 119) will be suppressed by ratios 
of a'j{k) / a"^, and thus the PBS perturbative expan- 
sion is well-behaved as long as we study scales larger 
than the Lagrangian size of the halos (which acts as 
the UV cutoff of the theory). In fact the resulting 
expansion is like in effective field theory, where the 
corrections are given by powers of the scale k over 
the cutoff RJ^^ ■ 

A different view on loops in PBS is presented 
in [29], where they argue that the PBS bias param- 



eters correspond to the renormalized bias parame- 
ters of [78, 85] and thus when computing loops one 
should discard terms which are UV sensitive. In 
this proposal, the running due to loops controlled 
by ('^^(fc)/'^™) is absent, so this makes distinct pre- 
dictions from our approach. Detailed calculations of 
correlation functions including loops and compari- 
son against simulations will be presented elsewhere. 



VI. CONCLUSIONS 

In this paper we have studied two major issues, 
the generation of initial conditions in N-body sim- 
ulations with primordial non-Gaussian (PNG) non- 
local models motivated by inflation and the peak- 
background split (PBS) method for predicting the 
scale-dependent bias in these models. The main re- 
sults in each case are, 

i) using as an example the factorizable tem- 
plates motivated by inflationary bispcctra, we 
showed how to construct initial conditions for 
nonlocal PNG models in an efficient way. Our 
algorithm is only slightly slower for nonlocal 
PNG than for the local case, and scales with 
particle number in the same way (ATpa^ In A^par) 
as Gaussian initial conditions. The result- 
ing implementation is given in Eqs. (34-35) 
for equilateral and orthogonal templates. Ap- 
pendix A generalizes these results beyond 
scale- invariance and Appendix B gives the gen- 
eral prescription for more general templates 
(including scale-dependent /nl models) and 
non-local ^nl PNG. 

ii) we derived predictions for the scale-dependent 
bias from the PBS in the context of the 
excursion-set approach to halos for general, 
non-local, PNG to describe scale-dependent 
bias. Our main result for the scale-dependent 
bias contribution in generic /nl PNG is given 
by Eq. (76), which generalizes previous results 
in the literature to go beyond the assiimption 
of universality of the mass function, Markovian 
evolution of the random walk, and includes be- 
yond leading-order effects due to the transfer 
function, filter effects and gradients. Our re- 
sults, which only require an accurate knowl- 
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edge of the Gaussian mass function, show im- 
proved agreement against numerical simula- 
tions when compared to the predictions pre- 
viously available in the literature (see Figs. 5- 
8). We also present, for the first time, results 
for the quadratic bias parameters for non-local 
PNG. 

Let us now expand on each of these in more detail. 

Regarding the generation of initial conditions, we 
showed that there arc multiple solutions for the non- 
local kernel given a desired bispectrum, but our gen- 
eral method does lead to a (sum of) factorizable ker- 
nel for a given factorizable bispectrum, preserving 
this advantage for a fast numerical generation of ini- 
tial conditions. We discussed in detail how to re- 
strict the linear combination of solutions using the 
one-loop corrections to the power spectrum to pre- 
serve the large-scale limit given by primordial corre- 
lations. Even after these restrictions, there is still a 
one-parameter family of solutions for the equilateral 
and orthogonal templates. This lack of uniqueness 
may even show up in the squeezed limit (see Ap- 
pendix E). This freedom is actually welcomed, as it 
should help reproduce the correct "snake topology" 
(proportional to tnl) of the trispectrum in these 
models, something that is left for future work. 

On the other hand, we stressed that while the free- 
dom may be removed by choosing the reduced bis- 
pectrum as a kernel [33] , doing so imposes a restric- 
tive symmetry (under permutations of k\^k2,kz) 
that is not clearly required, it imposes a non-trivial 
constraint between the power spectrum, bispectrum 
and snake amplitude of the four-point function that 
is unlikely to hold, and making the kernel non- 
factorizable complicates substantially the generation 
of initial conditions. 

We then used our algorithm to generate initial 
conditions to run a large suite of numerical sim- 
ulations (12 realizations of 1280^ particles in a 
2Ah~^ Gpc box) with PNG initial conditions of lo- 
cal, equilateral and orthogonal type and study the 
scale-dependent bias of dark matter halos. 

Concerning the implementation of PBS for generic 
PNG, we find that our predictions for the amplitude 
of the scale dependence are in excellent agreement 
with the measurements in our numerical simulations 
for a variety of halo masses, redshifts and halo def- 
initions. Our prediction for /nl PNG, Eq. (76), is 



valid for any non-local PNG and is simple to imple- 
ment (requiring only knowing the mass function for 
Gaussian simulations with the same cosmology) and 
we demonstrated it to be more accurate than previ- 
ous results in the literature that assume universal- 
ity (Eq. 39), Markovian evolution (Eq. 44) and the 
large-scale limit {k 0). Therefore, with our im- 
proved predictions, there is no need to include any 
fudge factors to describe scale-dependent bias, and 
the result can be written for /nl PNG in terms of 
the Gaussian mass function independent of the dy- 
namics assumed for halo formation. In addition, we 
extended our PBS formalism to deal with non-local 
5NL PNG. 

We also stressed that in doing the PBS the choice 
of the split variable (the Gaussian vs the non- 
Gaussian $) may lead to difi^erent theoretical predic- 
tions. For local and orthogonal models constraining 
the large-scale modes using or gives consistent 
answers at ©(/nl) at large scales (though they lead 
to sub-leading differences at small scales described 
by the difference between the cross bispectrum B^^^ 
and 5$$$). However, under special circumstances 
(when the long-mode kernel is sufficiently singular) 
these two predictions may differ even in their scal- 
ings in the squeezed limit (see Appendix E) and lead 
to different predictions for the scale-dependent bias 
at 0(/nl)- However this is not surprising, since in 
this case cannot be considered a function of just 
(j)e, as (/)s-modes can give an important contribution 
to ^e, and thus constraining the large-scale (p and $ 
can differ significantly. 

An example of this situation is given by our implc^- 
mentation of the equilateral template kernel, where 
a k~^^ bias is predicted in the PBS in ^ while a scale- 
independent bias is predicted by doing the split in 
In practice, the difference in the predictions is 
however very small (2% in the most biased samples) 
and thus unlikely to be important in practice and 
inconclusive on which approach is the preferred one. 
Wo pointed out, in addition, that the choice of the 
split variable can however be important for the halo 
bispectrum even in models where the kernels are reg- 
ular (as in local PNG), as imposing the large-scale 
constraint in (j> vs the non-Gaussian $ leads to large- 
scale different predictions for quadratic bias param- 
eters at 0{f^jj. This can be probed by measuring 
the halo bispectrum as a function of triangle shape, 
which we will address elsewhere. 
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In addition, in Section V wc showed that the pre- 
diction of perturbative local bias models that the 
amplitude of the scale-dependent bias is controlled 
by the quadratic bias 62 does not agree with simula- 
tions, where halos with negative 62 have an enhance- 
ment of bias at low-A: rather than suppression. 

Finally, we considered what happens with loop 
corrections due to nonlinear bias in the PBS, and 
showed that these are suppressed by the large- to 
small-scale variance aj{k)/a^, and thus they in- 
duced runnings with k away from the tree-level 
predictions that are well under control as long as 
fci?L <C 1, where i?L is the Lagrangian size of ha- 
los of mass m. As this scale is approached, bias 
ceases to be local (in UV sense) and the perturba- 
tive expansion breaks down due to the effects of halo 
exclusion. However, this expansion should be well- 
behaved when studying halos and galaxies at large- 
scalcs, which is where most of the constraining power 
on PNG lies. We will report on this in detail in the 
near future. 

While this work was being prepared for submis- 
sion, the paper [88] appeared in which the algorithm 
we present here (in the generic form given in Ap- 
pendix B) is also put forward to construct non-local 
PNG fields. They, however, used a different method 
to generate the kernel (the symmetrized kernel given 
by the reduced bispectrum). 



The parallel code using 2LPT initial conditions 
for non-local PNG models developed for this work is 
available^. 
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Appendix A: Generating Non-Local Primordial 
non-Gaussianity for non-scale invariant models 

When the power spectrum is not scale-invariant, 
the non-local operators used in the main text, 



where m = ±1, ±2/3, ±1/3 for the templates we 
have used in this paper. The kernels in real space 
then read. 



ireq[0, <l>] = -3(1 - U) 02 - 2p2/3([^-l/3'/']') ± (4 " 3m) Pi/3(0[P_i/30]) ± 2 P2/3(0[P-2/3<^]) 

-3uVi{(f>[V-i4'] - ['P-2/3<t>] [P-i/3<f>]) (A2) 

and 

ifort[<^, 0] = -9(1 -U)<l)^- 8p2/3([^-l/3</']') ± (10 - 9u) Pi/3((/'[P-l/3<^]) + 8 ^2/3 ('/'[7'-2/3'^]) 

-9uVi{<f>[V-i4'] - [P-2/3<t>] [P-i/s<l>]) (A3) 

For scale-invariant power spectra, ^±2/3 = — V^^ and ^±1/3 = d^^ (Eqs. 13-15), which yields the results 
quoted in the text (Eqs. 34-35). As discussed in the text, wc set the free parameter u to zero for simplicity. 
To minimize sub-leading one-loop corrections in the power spectrum that go as k~^ one may instead choose 
u as given in Eq. (33). Other possible choices are discussed in Appendix E. 

I 



Appendix B: More General Initial Conditions, 
Including Scale-Dependent /nl and 5nl 

In this Appendix, we give a short, self-contained 
recipe for generating initial conditions for fairly gen- 
eral, but factorizable, primordial non-gaussianity. 
Let's divide the discussion into two parts. One is 
how to simulate efhciently a kernel of the factoriz- 
able type. The other is how to determine the right 
kernel to use for a bispectrum of the factorizable 
type. 

First, how to simulate. Suppose one wants to sim- 



ulate: 

$(x) = (Pix) + j d'hiffka K12 0(ki)0(k2) 

with a kernel K12 of the form: 

K,2 = J2D /i(fci) /2(fc2) /3(|ki + k2l) (B2) 

where /i, /2 and /a are arbitrary functions, and the 
summation ^ denotes the fact that the kernel could 
be the sum of many terms of such factorizable form 
with different functions /i, /2, /s and constants D. 
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Our algorithm woiild allow us to simulate each of 
these, and one can simply add to obtain the desired 
non-Gaussian Note that we have subsumed /nl 
into the definition K12 to allow for the fact that even 
/nl could be scale dependent. 

The simulation algorithm goes as follows. 
First, generate a Gaussian random in Fourier 
space. Then, Fourier transform 4>(i<-i)fi{ki) and 
^(k2)/2(^2) separately back into real space and mul- 
tiply the results: let's denote this product in real 
space by [/i(^](x) x [f2(j)](x). Then, Fourier trans- 
form this back to Fourier space and multiply by /a, 
and finally Fourier transform back: 

D J d'he-'^^^-^Mh) 

I d^ye^'^^-yiMKy) x [M]{y) (B3) 

This real space quantity successfully realizes 
J d^kid^k2Ki24>{ki)4){k2)e-^^^^+^^>'' with the ker- 
nel given by K12 = D/i(fci)/2(fc2)/3(|ki + kal). 
Note that the algorithm requires only Fourier trans- 
forms and multiplications (in Fourier or real space). 
Nowhere do we need to perform expensive convolu- 
tions. 

Now, for a given bispectrum of the factorizable 
type, how do we decide what kernel to use? Suppose 
one wants to simulate a bispectrum of the form 

-B123 = '^^C[gi{ki)g2{k2)g3{k3) + 5 perm.] 

(B4) 

where gi , 32 and 53 are arbitrary functions, and the 
summation J2 denotes the fact that one might have 
a bispectrum that is a sum over terms of the above 
form, but with different functions gi, .92, .93 and con- 
stants C ^. Once again, if we know how to deter- 
mine the kernel for one of these, we can simply add 
to obtain the total. A kernel that produces the bis- 
pectrum 2C[gi{ki)g2{k2)g3{k3) + 5 perm.] is: 

K^2 = C^-^ [gi{ki)g2{k2)+g2{ki)gM 

(B5) 



Even non-separable bispectra can actually be well approx- 
imated in this way by using an appropriate basis [37, 39]. 



where Pi and P2 are the (p power spectrum at ki 
and k2- The kernel is thus of the form in Eq. (B2), 
and can therefore be simulated using the FFT-based 
algorithm above. 

This is not the only choice of kernel that will pro- 
duce the desired bispectrum. There are in fact two 
other kernels possible, which involve permutations 
of 1, 2 and 3, and the most general kernel is a linear 
combination of the three. However, for simplicity, 
one could adopt the kernel described above, but take 
care to check that the resulting one-loop <I> power 
spectrum will not have an IR limit that is more 
divergent than the tree level one (see Section II D 
and Appendix E for examples). This corresponds to 
choosing a g^ (since there are 3 functions involved, 
one has the freedom to choose which to call gs) such 
that (?3(fc)^/P^(A;) is not IR divergent. 

Finally, we note that following the same algorithm 
it is straightforward to construct non-local PNG ini- 
tial conditions for ^nl models. That is, a term in 
the Bardeen potential at cubic order in Eq. (Bl) of 
the form 

J d^krd^k2 L123 ^(ki)0(k2).^(k3)e-'(''i+''-+'^«)"' 

(B6) 

where L is the cubic kernel (proportional to ^nl) 
and assumed to be a sum of factorizable terms 

il23 = ^£^/l(fcl)/2(fc2)/3(fc3)/4(|ki +k2 +k3|) 

(B7) 

can also be easily implemented following the pro- 
cedure outlined above for non-local /nl-PNG. For 
work on local ^nl N-body simulations see [89, 90] . 

Appendix C: Further Tests on Initial Conditions 

We have done a number of tests on our initial con- 
ditions generator. The most direct test is to make 
sure that the bispectrum of the Bardeen potential is 
indeed the one desired for a given template. Here 
we do this for each of the operators that enter into 
the local, equilateral and orthogonal model. Let's 
consider the local model, for which we used 

$A = '^ + W (CI) 

and the corresponding bispectrum 

Ba = VMP^P2 + P2P3 + P3P1) (C2) 
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$D = + /nlV-2(^V2,/,) (C7) 
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FIG. 10. The ratio between the Bardeen potential Bis- 
pectrum measured in the initial conditions and the tree- 
level prediction for all triangles with sides smaller than 
60 times the fundamental mode for the four different 
operators discussed in the text (with amphtudes set by 
/nl ~ —400), from top to bottom: A (black), B (red), 
C (green) and D (blue), respectively. All ratios are close 
to one but are shifted by integer units for clarity. 

the second generator discussed in Section II B 

$B = <^ + /NLV-2(a0)2 (C3) 

for which the bispectrum reads 

Bb = -6/nl(PiP2P3)'/' (C4) 

We also consider the two first operators correspond- 
ing to the third generator in Section II C (see Eq. 26), 
i.e. 

*c = + /NLa-H05</') (C5) 

with bispectrum 

Be = /nl(Pi^'P2^'^3 + eye.) (C6) 



bispectrm 

Bd = ./nl(Pi /'P2^'-P3 + eye.) (C8) 



TABLE II. Best ht values of /nl obtained by htting the 
generated _B# to tree-level predictions for each of the four 
operators. 



Input /nl 


A 


B 


C 


D 


-50 


-48.5 


-48.0 


-48.4 


-48.4 


-400 


-390 


-397 


-392 


-392 


-3000 


-3224 


-3042 


-3078 


-3087 



The last operator in Eq. (26) is not included in these 
tests as we choose u = 2s = 0, and thus its ampli- 
tude vanishes (i.e. it does not appear in Eqs. 34 
or 35). 

To beat down cosmic variance we generated 100 
realizations each of initial conditions for these four 
operators with /nl — —50, —400, —3000 and mea- 
sured the bispectrum for each of them. A simple 
fit of the ratio of the numerically generated to 
its tree-level predictions (see Fig. 10, corresponding 
to the /nl — —400 case), gave the values shown in 
Table II for /nl, when using all triangles with sides 
up to fc 0.16 liMpc"^ corresponding to sixty times 
the fundamental mode of the box. The best fit val- 
ues are close to the input values even for large values 
of /nl where loop corrections may start to become 
important. 



Appendix D: New scale-independent corrections for local PNG 



While it is often assumed that for local PNG of /nl type the change in the small-scale density perturbation 
Ss corresponds to a local rescaling of the Gaussian small-scale matter fluctuations, Sf^ = (1 + 2/nl'7^«) ^f, 
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this is only true asymptotically at large-scales. In reality there is no local rescaling neither in real nor Fourier 
space, because of three different effects encoded by I21, which give quadratic corrections in {ke/kg) due to, 

1. contributions from the second term in V^($ — (f)) = 2/nl(?!'V^^ + | V<?i)p), which in previous literature 
has been neglected on the (incorrect) assumption that = for halos, 

2. contributions from the transfer function T{\k — 

3. contributions from the smoothing kernel Wrndfe — q\)- 

All of these are included in the main text exactly, by computing I12 numerically. Here we separate each 
of these contributions analj^ically by doing a Taylor expansion. The transfer function expansion in powers 

of {k/q) reads, 

T{\k-q\)= T{q) x 1 1 - (|-|) <5n + (^) ' [Sn [1 - 2(fe • qf] + (fe • qf{8n + 5nV2)] | (Dl) 

where 5n and 5n are the contributions to the effective spectral index and its running coming from the transfer 
function, 

5n{q) = Uefiiq) - Sn{q) = (D2) 

where for simplicity we have assumed a constant rig (no primordial running), and as usual HeS = dlnP/dlnk, 
with P the linear density power spectrum. We can use this to calculate, for example, I21 as a large-scale 
expansion in powers of k. For local PNG we have that the angular integral over k - q is given by 

/ (S)^™^''' " ^'^^ ^'^ ~ ^^'^ " '^'^ " '^^ WTK{qR)WMkR) T{q) 

<1 + 3"^ + ^ + 3"«^1 + + 1^ + + I2J) 

(D3) 



where Wk = d\nWTH{kR)/d\nk and we have used the standard summation theorem of Bessel functions to 
integrate over the top-hat window function (sec Appendix C in [35]). The cross-term k ■ q inside (fc — qf 
takes into account the |V^p contributions mentioned above. Since we are interested in small-fc limit we can 
use that Wk ~ —k'^R'^/5 and WTuikR) ~ 1 + Wfc/2 and thus, 

l2iik,m)^4f^^xi^al(l+lw,)+k' J ^P,(,)[M„(g)]^ [l + ^«,,(l+^n/2) + ^ + i^ + g] | (D4) 

We see from this that to zero-th order in k, we recover Eq. (81). Since the scale-dependent bias goes as 
I2i/M{k) oc (see Eq. 76), the next-to-leading order terms contribute to scale-independent bias, with 

a coefficient determined by integrating the velocity power spectrum against contributions from filtering and 
transfer function effects. 

I 

Appendix E: Squeezed Limit of Kernels and variable should be the Gaussian field (j) (as it has 
Scale-Dependent Bias been done in the past), or $ (as we discuss in Sec- 



Recall that in the main text we raised the issue 
of whether when implementing the PBS the split 



tion IIIF), which may give rise to different scai- 
ings even in the squeezed limit (agreeing with local 
bias models). These differences may arise because 
the former depends on the cross-bispectrum (-B$$<^) 
while the latter on the full bispcctrum (_B$$(i,), but 
how different these are from each other depends on 
the details of the long-mode kernel in the squeezed 
limit. In this appendix we want to explore these 
subtleties. 

Let's look at the squeezed limit of the kernels 
(which determine the low-fc limit) for equilateral and 
orthogonal PNG, Eqs. (29-30), before and after im- 
posing one-loop correction constraints. First, equi- 
lateral. There are actually two kernels that enter, de- 
pending on whether the coupling between fci and k2 
is to a large-scale {ke) or small-scale mode {kg '> kg), 

<Hfe.,fe,)«^(«-2s) + (l-3t)(|) 

+3(s + t)(|)%... (El) 

and 

KW{ks,K) « -3(«- 2s)(|)' - 2(1 - 3i)(|)' 

-6(. + t-l)(|)+... 

+ ... (E2) 

where \ks + k'J = k^. In the first kernel, to sim- 
plify the expressions we assumed that the third side 
is ks (i.e. Eq. El is for a fixed angle between kg 
and ki, the full expression as a function of angle is 



given below after one-loop correction constraints are 
imposed). Equation (E2), on the other hand, is for 
\k'g\ = kg and third side equal to kc, in this case the 
angle between ks and k'^ is restricted to be nearly 
TT to make a low-fc mode. Thus we are considering 
precisely the same triangle for both kernels, to high- 
light that they can differ significantly, even in the 
squeezed limit. 

We see from these two expressions that these 
kernels can have very different behavior, and that 
assuming that kernels arc totally symmetric in 
ki,k2, is very restrictive, since it forces equality of 
these two expressions. Even after imposing one-loop 
constraints the behavior of these kernels will gener- 
ically differ from each other in the squeezed limit. 

Before we get into one-loop corrections, note that 
requiring that ($(fe)) = leads to 

fc(fc) lim f Ki^^{ks,K)P4k,)d^ks (E3) 

ke^O J 

which from Eq. (E2) will typically diverge in the 
infrared (IR). However, this IR divergence is no dif- 
ferent that in the local case, in which Eq. (E3) leads 
to the variance of the </> field, which is also IR di- 
vergent, and can be cured by putting the field in 
a box, i.e. imposing that Pcj,{k) vanish for k < 
e = 27r/_Lbox- The same is true here, because the 
factors of (ks/ke)" in Eq. (E2) actually arise from 
P<i>{k() / P^{ks) ratios when solving for the kernels 
and thus cutting the P^{kf) power spectrum below 
some IR cutoff e regulates (<I>(fe)) in the same way 
as for local models. Therefore, we won't impose any 
constraint on the free parameters of the kernel from 
these considerations. 

From Eqs. (E1-E2) we can compute the squeezed 
limit of the relevant equilateral PNG bispectra, 



Bll^iks, fc„ ke) « /nl [2if (^)(fe,„ kt)P^{ke)P4k,) + 2K(^)(fc;, ke)P^{ke)P4k',)\ , (E4) 

which depends on the parameters s, t, u and 

Bl\^{k,,k,M « /nl [2Ki:^{KM)P^{kf:)P^{k,) + 2Ki:^{k'MP<t>ikt)PM) + 2^iqHfcs, K)[P'^i.k,) 

'ki^^ 



= 12/nl(^) P^{ki)P^{K)^... 



(E5) 



which, of course doesn't depend on s, t, u, as varying these parameters leaves the tree-level bispectrum of 
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$ unchanged. Note that the relevant bispectrum for 
the calculation of the bias when doing the PBS in 
(f) is -Bjj0 = [Mm{ks)]'^ B^^cf,, thus we are interested 
in the squeezed limit of -6$$^ and how it relates to 
that in We see from Eq. (El) that there is 

a priori no reason for Eqs. (E5) and (E4) to agree, 
unless the first two terms in Eq. (El) vanish and the 
last one has a suitable coefRcient. 

Let us now impose the constraints from one-loop 
corrections, which are easy to see from Eq. (E2) since 
that's the kernel whose square enters into the cal- 
culation of the one-loop power spectrum. We see, 
in agreement with the detailed calculations of the 



power spectrum presented in Section II D, that to 
avoid IR corrections more important than k^^ we 
must require u = 2s (to avoid kj^) and t = 1/3 (to 
avoid kj^), which then says that 

if«(*..M-f(»+|)(|)%... (E6) 

and 

ifW(fc„fc;)«(4-3u)(^)+... (E7) 

Let us now restore the full angular dependence in 
Eq. (E6), and include the cosine of the angle between 
short and long modes x = kg ■ ke, 



Ki'^\ks, ki) « 2(-2 + 3w) x(^) + [{7x^ + 2x-l)- '^u i9x^ + 2x- 3)] (^)' + . . . 



(E8) 



from which Eq. (E6) can be obtained by using x = —k(/2ks for the triangle considered above. We see from 

Eq. (E8) that a term of 0{ki/ks) is now present, it is however proportional to the cosine variable x, which 
as long the triangle is not isosceles is of order unity. This means that the squeezed limit in Eq. (E4) is given 

by 



/NL[4-12x2+6«(l+a;2)] (^)' x P^(A;^)P^(fc,) + 



(E9) 



where the linear term in {ki/kg) proportional to 
X cancels when summing over the two terms in 
Eq. (E4) since kg and fc^ are nearly antiparal- 
lel. On the other hand, restoring the full angular 
dependence for the Bardeen potential bispectrum, 
Eq. (E5), we have 

« 12/nl (1 - x^) (^^y P^{ke)P^{ks) 

(ElO) 

We see that Eqs. (E9) and (ElO) disagree in ampli- 



tude for all values of u, although they agree on the 
scaling. Does this mean that PBS predicts scale- 
independent bias corrections for the equilateral tem- 
plate? It's not so trivial, because Eq. (E9) depends 
on u, and in fact for our choice here {u = 0) the 
angular average of the (ke/kg)^ amplitude vanishes! 
Thus, there is no correction to the bias at this or- 
der, and the leading order will come from terms of 
order {ke/ks)^. The cross bispectrum to this order 
for u = reads. 



^$<i>0l«=o ~ /nl 



(4- 12a;^ 



^)(^)' + 2:r(21a;2 + 2ar-ll)( 



X P^iki)P^iks) + 



(Ell) 



The full calculation to this order will have to in- elude the cross term between the (fc^/fcs)^ and the 
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next-to-leadiiig corrections 0{k(,/ks) coming from 
Mm{\k — q\) discussed in Appendix D. However, 
ignoring such cross terms and using only the angu- 
lar average of the second term in Eq. (Ell) gives a 
very good approximation to the asymptotic value of 
/21 in the low-A; limit as shown by the bottom panel 
in Fig. 3. 



Therefore we have an interesting and nontrivial 
example for the equilateral template where the scal- 
ing of the kernel, Eq. (E8), would suggest a k~^ 
scale- dependent bias, this term however cancels in 
the cross bispectrum, Eq. (E9), which suggests a 
scale-independent bias, but this amplitude has for 
our choice u = a vanishing angular average, there- 
fore one has to go to the; next order and finds a 
scale- dependent bias contribution going as A;"*"^ with 
a small amplitude. Therefore, in this particular ex- 
ample, the scaling of the possibly scale-dependent 
bias when doing the PBS in (fc'*'^) disagrees from 
when doing the PBS in $ which agrees with 

local bias models. This is made possible because 
K^^^ is sufficiently singular compared to /C'^*' in the 
squeezed limit, i.e. K^^^ <^ K^^\ks/ke)^, to compen- 
sate for the ratio of power spectra between short and 
long modes. These differences could be in principle 
be checked in simulations by running initial condi- 
tions with different values of u. However, because 
for equilateral PNG these corrections are small, it 
may be difficult to do in practice (see Fig. 8). 



Let us now turn to the orthogonal template. We 



have for the kernels. 

Kin (ks, k,) « ^ (,, - 2s) + (1 - 90 (|) + . . . 

(E12) 

and 

i^W(fe„fe'J « -9(« - 2.)(|)' - 2(4 - 9t)(|)' 
-18(s + f-l)(^) +... (E13) 

which after imposing one-loop power spectrum con- 
straints (u = 2s and t = 4/9) give, restoring the 
angular dependence 

(E14) 

and 

K^J^lik,, fe'J « [10 - 9^.(1 + x^)] (|) + . .(E15) 

again, showing very different scaling of the kernels. 
However, because the K'-^^ kernel is only enhanced 
by (ks/ke)^ compared to K^^^ it cannot compensate 
for the suppression factor P^{ks) / Ptpiki) in the <I> 
bispectrum and thus now both squeezed-limit bis- 
pectra agree to leading order, 

BT^^ « BJ'J* « -12/NL(^)j'0(fc«)P4^.)' 

(E16) 

where again the linear terms in x in Eq. (E14) cancel 
in the cross bispectrum, leading to a u-independent 
result. In this case, therefore, the scaling of the 
scale-dependence bias with k^ would agree in the 
PBS with the split performed in (j) or 
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